From ac152578c237aa3796db7bff175f07005123404d Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Wed, 27 Oct 2010 20:45:07 +0800 Subject: [PATCH] clean up: more bug fixes and added intermediate, type specific factorization classes --- src/MathNet.Numerics.5.1.ReSharper | 6 +- .../Complex/Factorization/Cholesky.cs | 81 ++ .../Complex/Factorization/DenseCholesky.cs | 45 +- .../Complex/Factorization/DenseEvd.cs | 18 +- .../Complex/Factorization/DenseGramSchmidt.cs | 36 +- .../Complex/Factorization/DenseLU.cs | 25 +- .../Complex/Factorization/DenseQR.cs | 36 +- .../Complex/Factorization/DenseSvd.cs | 37 +- .../Complex/Factorization/Evd.cs | 114 +++ .../Complex/Factorization/GramSchmidt.cs | 89 ++ .../LinearAlgebra/Complex/Factorization/LU.cs | 68 ++ .../LinearAlgebra/Complex/Factorization/QR.cs | 91 ++ .../Complex/Factorization/Svd.cs | 119 +++ .../Complex/Factorization/UserCholesky.cs | 36 +- .../Complex/Factorization/UserEvd.cs | 18 +- .../Complex/Factorization/UserGramSchmidt.cs | 28 +- .../Complex/Factorization/UserLU.cs | 17 +- .../Complex/Factorization/UserQR.cs | 29 +- .../Complex/Factorization/UserSvd.cs | 34 +- .../Complex32/Factorization/Cholesky.cs | 81 ++ .../Complex32/Factorization/DenseCholesky.cs | 45 +- .../Complex32/Factorization/DenseEvd.cs | 18 +- .../Factorization/DenseGramSchmidt.cs | 36 +- .../Complex32/Factorization/DenseLU.cs | 25 +- .../Complex32/Factorization/DenseQR.cs | 36 +- .../Complex32/Factorization/DenseSvd.cs | 37 +- .../Complex32/Factorization/Evd.cs | 116 +++ .../Complex32/Factorization/GramSchmidt.cs | 89 ++ .../Complex32/Factorization/LU.cs | 68 ++ .../Complex32/Factorization/QR.cs | 91 ++ .../Complex32/Factorization/Svd.cs | 119 +++ .../Complex32/Factorization/UserCholesky.cs | 37 +- .../Complex32/Factorization/UserEvd.cs | 18 +- .../Factorization/UserGramSchmidt.cs | 28 +- .../Complex32/Factorization/UserLU.cs | 17 +- .../Complex32/Factorization/UserQR.cs | 29 +- .../Complex32/Factorization/UserSvd.cs | 34 +- .../Double/Factorization/Cholesky.cs | 81 ++ .../Double/Factorization/DenseCholesky.cs | 45 +- .../Double/Factorization/DenseEvd.cs | 18 +- .../Double/Factorization/DenseGramSchmidt.cs | 35 +- .../Double/Factorization/DenseLU.cs | 25 +- .../Double/Factorization/DenseQR.cs | 36 +- .../Double/Factorization/DenseSvd.cs | 37 +- .../LinearAlgebra/Double/Factorization/Evd.cs | 114 +++ .../Double/Factorization/GramSchmidt.cs | 88 ++ .../LinearAlgebra/Double/Factorization/LU.cs | 67 ++ .../LinearAlgebra/Double/Factorization/QR.cs | 90 ++ .../Double/Factorization/SparseCholesky.cs | 258 ----- .../Double/Factorization/SparseLU.cs | 317 ------ .../Double/Factorization/SparseQR.cs | 358 ------- .../Double/Factorization/SparseSvd.cs | 950 ------------------ .../LinearAlgebra/Double/Factorization/Svd.cs | 118 +++ .../Double/Factorization/UserCholesky.cs | 37 +- .../Double/Factorization/UserEvd.cs | 18 +- .../Double/Factorization/UserGramSchmidt.cs | 27 +- .../Double/Factorization/UserLU.cs | 18 +- .../Double/Factorization/UserQR.cs | 29 +- .../Double/Factorization/UserSvd.cs | 34 +- .../Generic/Factorization/Cholesky.cs | 91 +- .../Generic/Factorization/Evd.cs | 130 +-- .../Generic/Factorization/GramSchmidt.cs | 41 +- .../LinearAlgebra/Generic/Factorization/LU.cs | 113 +-- .../LinearAlgebra/Generic/Factorization/QR.cs | 93 +- .../Generic/Factorization/Svd.cs | 101 +- .../Single/Factorization/Cholesky.cs | 81 ++ .../Single/Factorization/DenseCholesky.cs | 45 +- .../Single/Factorization/DenseEvd.cs | 18 +- .../Single/Factorization/DenseGramSchmidt.cs | 36 +- .../Single/Factorization/DenseLU.cs | 25 +- .../Single/Factorization/DenseQR.cs | 36 +- .../Single/Factorization/DenseSvd.cs | 37 +- .../LinearAlgebra/Single/Factorization/Evd.cs | 115 +++ .../Single/Factorization/GramSchmidt.cs | 88 ++ .../LinearAlgebra/Single/Factorization/LU.cs | 67 ++ .../LinearAlgebra/Single/Factorization/QR.cs | 90 ++ .../LinearAlgebra/Single/Factorization/Svd.cs | 118 +++ .../Single/Factorization/UserCholesky.cs | 38 +- .../Single/Factorization/UserEvd.cs | 18 +- .../Single/Factorization/UserGramSchmidt.cs | 30 +- .../Single/Factorization/UserLU.cs | 17 +- .../Single/Factorization/UserQR.cs | 29 +- .../Single/Factorization/UserSvd.cs | 36 +- src/Numerics/Numerics.csproj | 28 +- src/Silverlight/Silverlight.csproj | 76 +- .../Complex/Factorization/EvdTests.cs | 22 +- .../Complex/Factorization/UserEvdTests.cs | 22 +- .../Complex32/Factorization/EvdTests.cs | 24 +- .../Complex32/Factorization/UserEvdTests.cs | 24 +- .../Complex32/MatrixTests.Arithmetic.cs | 36 +- .../Double/Factorization/EvdTests.cs | 22 +- .../Double/Factorization/UserEvdTests.cs | 22 +- .../Double/MatrixTests.Arithmetic.cs | 36 +- .../Single/Factorization/EvdTests.cs | 24 +- .../Single/Factorization/UserEvdTests.cs | 22 +- .../Single/MatrixTests.Arithmetic.cs | 36 +- .../LinearAlgebraTests/Single/MatrixTests.cs | 24 +- 97 files changed, 2719 insertions(+), 3843 deletions(-) create mode 100644 src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs create mode 100644 src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs create mode 100644 src/Numerics/LinearAlgebra/Complex/Factorization/GramSchmidt.cs create mode 100644 src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs create mode 100644 src/Numerics/LinearAlgebra/Complex/Factorization/QR.cs create mode 100644 src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs create mode 100644 src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs create mode 100644 src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs create mode 100644 src/Numerics/LinearAlgebra/Complex32/Factorization/GramSchmidt.cs create mode 100644 src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs create mode 100644 src/Numerics/LinearAlgebra/Complex32/Factorization/QR.cs create mode 100644 src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs create mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs create mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs create mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/GramSchmidt.cs create mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/LU.cs create mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/QR.cs delete mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/SparseCholesky.cs delete mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/SparseLU.cs delete mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/SparseQR.cs delete mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/SparseSvd.cs create mode 100644 src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs create mode 100644 src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs create mode 100644 src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs create mode 100644 src/Numerics/LinearAlgebra/Single/Factorization/GramSchmidt.cs create mode 100644 src/Numerics/LinearAlgebra/Single/Factorization/LU.cs create mode 100644 src/Numerics/LinearAlgebra/Single/Factorization/QR.cs create mode 100644 src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs diff --git a/src/MathNet.Numerics.5.1.ReSharper b/src/MathNet.Numerics.5.1.ReSharper index 778ad35b..3bcec989 100644 --- a/src/MathNet.Numerics.5.1.ReSharper +++ b/src/MathNet.Numerics.5.1.ReSharper @@ -25,7 +25,11 @@ indices Frobenius Pointwise multipcation -kronecker +kronecker +Cholesky +Eigen +mxn +nxn diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs new file mode 100644 index 00000000..435540d7 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs @@ -0,0 +1,81 @@ +// +// 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.Complex.Factorization +{ + using System.Numerics; + using Generic.Factorization; + + /// + /// A class which encapsulates the functionality of a Cholesky factorization. + /// For a symmetric, positive definite matrix A, the Cholesky factorization + /// is an lower triangular matrix L so that A = L*L'. + /// + /// + /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric + /// or positive definite, the constructor will throw an exception. + /// + public abstract class Cholesky : Cholesky + { + /// + /// Gets the determinant of the matrix for which the Cholesky matrix was computed. + /// + public override Complex Determinant + { + get + { + var det = Complex.One; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; + } + + return det; + } + } + + /// + /// Gets the log determinant of the matrix for which the Cholesky matrix was computed. + /// + public override Complex DeterminantLn + { + get + { + var det = Complex.Zero; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det += 2.0 * CholeskyFactor[j, j].NaturalLogarithm(); + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs index 69ccf087..8cf98d70 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; using Threading; @@ -46,7 +45,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 class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -111,13 +110,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -160,13 +159,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -176,39 +175,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dfactor = (DenseMatrix)CholeskyFactor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex AddT(Complex val1, Complex val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override Complex LogT(Complex val1) - { - return val1.NaturalLogarithm(); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs index 94820088..188a8d45 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,16 +47,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public DenseEvd(DenseMatrix matrix) { @@ -949,16 +948,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs index 048b84bc..eeee8398 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; using Threading; @@ -44,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class DenseGramSchmidt : GramSchmidt + public class DenseGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an unitary matrix @@ -154,13 +153,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -199,41 +198,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs index 058e4742..0233dd39 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; using Threading; @@ -45,7 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class DenseLU : LU + public class DenseLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -112,13 +111,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -161,13 +160,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -188,19 +187,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization Control.LinearAlgebraProvider.LUInverseFactored(result.Data, result.RowCount, Pivots); return result; } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs index 1fa36654..22dcd05f 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -45,7 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class DenseQR : QR + public class DenseQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -110,13 +109,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -155,41 +154,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs index aa98288e..783df203 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -49,7 +48,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,7 +56,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public DenseSvd(DenseMatrix matrix, bool computeVectors) { @@ -118,13 +117,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, input.ColumnCount, dresult.Data); @@ -168,40 +167,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs new file mode 100644 index 00000000..77c63a7a --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/Evd.cs @@ -0,0 +1,114 @@ +// +// 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.Complex.Factorization +{ + using System.Numerics; + using Generic.Factorization; + + /// + /// Eigenvalues and eigenvectors of a real matrix. + /// + /// + /// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is + /// diagonal and the eigenvector matrix V is orthogonal. + /// I.e. A = V*D*V' and V*VT=I. + /// If A is not symmetric, then the eigenvalue matrix D is block diagonal + /// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, + /// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The + /// columns of V represent the eigenvectors in the sense that A*V = V*D, + /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly + /// conditioned, or even singular, so the validity of the equation + /// A = V*D*Inverse(V) depends upon V.Condition(). + /// + public abstract class Evd : Evd + { + /// + /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. + /// + public override Complex Determinant + { + get + { + var det = Complex.One; + for (var i = 0; i < VectorEv.Count; i++) + { + det *= VectorEv[i]; + + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + return 0; + } + } + + return det.Magnitude; + } + } + + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + var rank = 0; + for (var i = 0; i < VectorEv.Count; i++) + { + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + continue; + } + + rank++; + } + + return rank; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < VectorEv.Count; i++) + { + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/GramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/GramSchmidt.cs new file mode 100644 index 00000000..ffcff930 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/GramSchmidt.cs @@ -0,0 +1,89 @@ +// +// 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.Complex.Factorization +{ + using System; + using System.Numerics; + using Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization. + /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix. + /// + /// + /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. + /// + public abstract class GramSchmidt : GramSchmidt + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override Complex Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = Complex.One; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) + { + return 0; + } + } + + return det.Magnitude; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs new file mode 100644 index 00000000..47333fa0 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs @@ -0,0 +1,68 @@ +// +// 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.Complex.Factorization +{ + using System.Numerics; + using Generic.Factorization; + + /// + /// 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 : LU + { + /// + /// Gets the determinant of the matrix for which the LU factorization was computed. + /// + public override Complex Determinant + { + get + { + var det = Complex.One; + for (var j = 0; j < Factors.RowCount; j++) + { + if (Pivots[j] != j) + { + det *= -Factors.At(j, j); + } + else + { + det *= Factors.At(j, j); + } + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/QR.cs new file mode 100644 index 00000000..3c5d9161 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/QR.cs @@ -0,0 +1,91 @@ +// +// 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.Complex.Factorization +{ + using System; + using System.Numerics; + using Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition. + /// Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m) + /// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix + /// (also called right triangular matrix). + /// + /// + /// The computation of the QR decomposition is done at construction time by Householder transformation. + /// + public abstract class QR : QR + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override Complex Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = Complex.One; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) + { + return 0; + } + } + + return det.Magnitude; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs new file mode 100644 index 00000000..f321d76c --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs @@ -0,0 +1,119 @@ +// +// 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.Complex.Factorization +{ + using System; + using System.Linq; + using System.Numerics; + using Generic; + using Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the singular value decomposition (SVD). + /// Suppose M is an m-by-n matrix whose entries are real numbers. + /// Then there exists a factorization of the form M = UΣVT where: + /// - U is an m-by-m unitary matrix; + /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; + /// - VT denotes transpose of V, an n-by-n unitary matrix; + /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal + /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined + /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M. + /// + /// + /// The computation of the singular value decomposition is done at construction time. + /// + public abstract class Svd : Svd + { + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + return VectorS.Count(t => !t.Magnitude.AlmostEqual(0.0)); + } + } + + /// + /// Gets the two norm of the . + /// + /// The 2-norm of the . + public override Complex Norm2 + { + get + { + return VectorS[0].Magnitude; + } + } + + /// + /// Gets the condition number max(S) / min(S) + /// + /// The condition number. + public override Complex ConditionNumber + { + get + { + var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; + return VectorS[0].Magnitude / VectorS[tmp].Magnitude; + } + } + + /// + /// Gets the determinant of the square matrix for which the SVD was computed. + /// + public override Complex Determinant + { + get + { + if (MatrixU.RowCount != MatrixVT.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = Complex.One; + foreach (var value in VectorS) + { + det *= value; + if (value.Magnitude.AlmostEqual(0.0)) + { + return 0; + } + } + + return det.Magnitude; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs index 1b40a88a..56ea890f 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs @@ -45,7 +45,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 class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -221,39 +221,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization result[i] = sum / CholeskyFactor.At(i, i); } } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex AddT(Complex val1, Complex val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override Complex LogT(Complex val1) - { - return val1.NaturalLogarithm(); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs index df1302ba..3d73f2e1 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,16 +47,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public UserEvd(Matrix matrix) { @@ -944,16 +943,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserGramSchmidt.cs index 2ce704b4..ea8bb7c6 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserGramSchmidt.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class UserGramSchmidt : GramSchmidt + public class UserGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an unitary matrix @@ -250,30 +249,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs index 03e17611..17118574 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -44,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class UserLU : LU + public class UserLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -299,19 +298,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization return Solve(inverse); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs index dd64b30f..e4c9a382 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs @@ -34,7 +34,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System.Linq; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -46,7 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class UserQR : QR + public class UserQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -92,7 +91,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Generate column from initial matrix to work array /// /// Initial matrix - /// The firts row + /// The first row /// The last row /// Column index /// Generated vector @@ -329,29 +328,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs index 172de227..49556ca9 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -49,7 +48,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,7 +56,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public UserSvd(Matrix matrix, bool computeVectors) { @@ -718,14 +717,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization db = z; } - /// dded + /// /// Calculate Norm 2 of the column in matrix starting from row /// /// Source matrix /// The number of rows in /// Column index /// Start row index - /// Norm2 (Euclidean norm) of trhe column + /// Norm2 (Euclidean norm) of the column private static double Cnrm2Column(Matrix a, int rowCount, int column, int rowStart) { var s = 0.0; @@ -936,30 +935,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization result[j] = value; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs new file mode 100644 index 00000000..eef56c06 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs @@ -0,0 +1,81 @@ +// +// 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.Complex32.Factorization +{ + using Generic.Factorization; + using Numerics; + + /// + /// A class which encapsulates the functionality of a Cholesky factorization. + /// For a symmetric, positive definite matrix A, the Cholesky factorization + /// is an lower triangular matrix L so that A = L*L'. + /// + /// + /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric + /// or positive definite, the constructor will throw an exception. + /// + public abstract class Cholesky : Cholesky + { + /// + /// Gets the determinant of the matrix for which the Cholesky matrix was computed. + /// + public override Complex32 Determinant + { + get + { + var det = Complex32.One; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; + } + + return det; + } + } + + /// + /// Gets the log determinant of the matrix for which the Cholesky matrix was computed. + /// + public override Complex32 DeterminantLn + { + get + { + var det = Complex32.Zero; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det += 2.0f * CholeskyFactor[j, j].NaturalLogarithm(); + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs index 1525112b..204b1dc8 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; using Threading; @@ -46,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 DenseCholesky : Cholesky + public class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -111,13 +110,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -160,13 +159,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -176,39 +175,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dfactor = (DenseMatrix)CholeskyFactor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex32 AddT(Complex32 val1, Complex32 val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override Complex32 LogT(Complex32 val1) - { - return val1.NaturalLogarithm(); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs index dc40cac2..265aba55 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -49,16 +48,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public DenseEvd(DenseMatrix matrix) { @@ -953,16 +952,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs index 42d21502..5f86ddbb 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; using Threading; @@ -44,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class DenseGramSchmidt : GramSchmidt + public class DenseGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an unitary matrix @@ -154,13 +153,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -199,41 +198,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs index c5eefeae..61ad70ed 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; using Threading; @@ -45,7 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class DenseLU : LU + public class DenseLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -112,13 +111,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -161,13 +160,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -188,19 +187,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization Control.LinearAlgebraProvider.LUInverseFactored(result.Data, result.RowCount, Pivots); return result; } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs index 9f66ace9..64ef6526 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -45,7 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class DenseQR : QR + public class DenseQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -110,13 +109,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -155,41 +154,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs index 143d6aab..c7b5cb05 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -49,7 +48,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,7 +56,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public DenseSvd(DenseMatrix matrix, bool computeVectors) { @@ -118,13 +117,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, input.ColumnCount, dresult.Data); @@ -168,40 +167,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs new file mode 100644 index 00000000..b268eaaf --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/Evd.cs @@ -0,0 +1,116 @@ +// +// 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.Complex32.Factorization +{ + using System; + using System.Numerics; + using Generic.Factorization; + using Numerics; + + /// + /// Eigenvalues and eigenvectors of a real matrix. + /// + /// + /// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is + /// diagonal and the eigenvector matrix V is orthogonal. + /// I.e. A = V*D*V' and V*VT=I. + /// If A is not symmetric, then the eigenvalue matrix D is block diagonal + /// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, + /// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The + /// columns of V represent the eigenvectors in the sense that A*V = V*D, + /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly + /// conditioned, or even singular, so the validity of the equation + /// A = V*D*Inverse(V) depends upon V.Condition(). + /// + public abstract class Evd : Evd + { + /// + /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. + /// + public override Complex32 Determinant + { + get + { + var det = Complex.One; + for (var i = 0; i < VectorEv.Count; i++) + { + det *= VectorEv[i]; + + if (((Complex32)VectorEv[i]).AlmostEqual(Complex32.Zero)) + { + return 0; + } + } + + return new Complex32(Convert.ToSingle(det.Magnitude), 0.0f); + } + } + + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + var rank = 0; + for (var i = 0; i < VectorEv.Count; i++) + { + if (((Complex32)VectorEv[i]).AlmostEqual(Complex32.Zero)) + { + continue; + } + + rank++; + } + + return rank; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < VectorEv.Count; i++) + { + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/GramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/GramSchmidt.cs new file mode 100644 index 00000000..d183d7d8 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/GramSchmidt.cs @@ -0,0 +1,89 @@ +// +// 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.Complex32.Factorization +{ + using System; + using Generic.Factorization; + using Numerics; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization. + /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix. + /// + /// + /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. + /// + public abstract class GramSchmidt : GramSchmidt + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override Complex32 Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = Complex32.One; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) + { + return 0; + } + } + + return det.Magnitude; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs new file mode 100644 index 00000000..094b9f8d --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs @@ -0,0 +1,68 @@ +// +// 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.Complex32.Factorization +{ + using Generic.Factorization; + using Numerics; + + /// + /// 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 : LU + { + /// + /// Gets the determinant of the matrix for which the LU factorization was computed. + /// + public override Complex32 Determinant + { + get + { + var det = Complex32.One; + for (var j = 0; j < Factors.RowCount; j++) + { + if (Pivots[j] != j) + { + det *= -Factors.At(j, j); + } + else + { + det *= Factors.At(j, j); + } + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/QR.cs new file mode 100644 index 00000000..c899365e --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/QR.cs @@ -0,0 +1,91 @@ +// +// 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.Complex32.Factorization +{ + using System; + using Generic.Factorization; + using Numerics; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition. + /// Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m) + /// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix + /// (also called right triangular matrix). + /// + /// + /// The computation of the QR decomposition is done at construction time by Householder transformation. + /// + public abstract class QR : QR + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override Complex32 Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = Complex32.One; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) + { + return 0; + } + } + + return det.Magnitude; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs new file mode 100644 index 00000000..5e3bd0ec --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs @@ -0,0 +1,119 @@ +// +// 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.Complex32.Factorization +{ + using System; + using System.Linq; + using Generic; + using Generic.Factorization; + using Numerics; + using Properties; + + /// + /// A class which encapsulates the functionality of the singular value decomposition (SVD). + /// Suppose M is an m-by-n matrix whose entries are real numbers. + /// Then there exists a factorization of the form M = UΣVT where: + /// - U is an m-by-m unitary matrix; + /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; + /// - VT denotes transpose of V, an n-by-n unitary matrix; + /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal + /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined + /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M. + /// + /// + /// The computation of the singular value decomposition is done at construction time. + /// + public abstract class Svd : Svd + { + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + return VectorS.Count(t => !t.Magnitude.AlmostEqual(0.0f)); + } + } + + /// + /// Gets the two norm of the . + /// + /// The 2-norm of the . + public override Complex32 Norm2 + { + get + { + return VectorS[0].Magnitude; + } + } + + /// + /// Gets the condition number max(S) / min(S) + /// + /// The condition number. + public override Complex32 ConditionNumber + { + get + { + var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; + return VectorS[0].Magnitude / VectorS[tmp].Magnitude; + } + } + + /// + /// Gets the determinant of the square matrix for which the SVD was computed. + /// + public override Complex32 Determinant + { + get + { + if (MatrixU.RowCount != MatrixVT.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = Complex32.One; + foreach (var value in VectorS) + { + det *= value; + if (value.Magnitude.AlmostEqual(0.0f)) + { + return 0; + } + } + + return det.Magnitude; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs index a33a852d..f7414b41 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -45,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 UserCholesky : Cholesky + public class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -221,39 +220,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization result[i] = sum / CholeskyFactor.At(i, i); } } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex32 AddT(Complex32 val1, Complex32 val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override Complex32 LogT(Complex32 val1) - { - return val1.NaturalLogarithm(); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs index a86371b9..803f0a86 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -49,16 +48,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public UserEvd(Matrix matrix) { @@ -948,16 +947,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs index 662a0c7b..3973c6dc 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class UserGramSchmidt : GramSchmidt + public class UserGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an unitary matrix @@ -250,30 +249,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs index 7d805ced..864c94e3 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -44,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class UserLU : LU + public class UserLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -299,19 +298,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization return Solve(inverse); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs index e1c02048..9d8b21ea 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization using System; using System.Linq; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -46,7 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class UserQR : QR + public class UserQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -92,7 +91,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Generate column from initial matrix to work array /// /// Initial matrix - /// The firts row + /// The first row /// The last row /// Column index /// Generated vector @@ -329,29 +328,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.cs index 1c906540..d0c21b42 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserSvd.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -49,7 +48,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -57,7 +56,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public UserSvd(Matrix matrix, bool computeVectors) { @@ -718,14 +717,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization db = z; } - /// dded + /// /// Calculate Norm 2 of the column in matrix starting from row /// /// Source matrix /// The number of rows in /// Column index /// Start row index - /// Norm2 (Euclidean norm) of trhe column + /// Norm2 (Euclidean norm) of the column private static float Cnrm2Column(Matrix a, int rowCount, int column, int rowStart) { var s = 0.0f; @@ -936,30 +935,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization result[j] = value; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs new file mode 100644 index 00000000..871a2fad --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs @@ -0,0 +1,81 @@ +// +// 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 Generic.Factorization; + + /// + /// A class which encapsulates the functionality of a Cholesky factorization. + /// For a symmetric, positive definite matrix A, the Cholesky factorization + /// is an lower triangular matrix L so that A = L*L'. + /// + /// + /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric + /// or positive definite, the constructor will throw an exception. + /// + public abstract class Cholesky : Cholesky + { + /// + /// Gets the determinant of the matrix for which the Cholesky matrix was computed. + /// + public override double Determinant + { + get + { + var det = 1.0; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; + } + + return det; + } + } + + /// + /// Gets the log determinant of the matrix for which the Cholesky matrix was computed. + /// + public override double DeterminantLn + { + get + { + var det = 0.0; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det += 2 * Math.Log(CholeskyFactor[j, j]); + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index 32a60253..ef80191b 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -44,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 DenseCholesky : Cholesky + public class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -109,13 +108,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -158,13 +157,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -174,39 +173,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dfactor = (DenseMatrix)CholeskyFactor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override double AddT(double val1, double val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override double LogT(double val1) - { - return Math.Log(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs index 5ec5ea22..6684d923 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,16 +47,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public DenseEvd(DenseMatrix matrix) { @@ -1222,16 +1221,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs index 3b190b9f..e749680f 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; using Threading; @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class DenseGramSchmidt : GramSchmidt + public class DenseGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an orthogonal matrix @@ -153,13 +152,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -198,40 +197,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs index b0bcacbe..8b33dc87 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class DenseLU : LU + public class DenseLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -110,13 +109,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -159,13 +158,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -186,19 +185,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Control.LinearAlgebraProvider.LUInverseFactored(result.Data, result.RowCount, Pivots); return result; } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs index 0d559497..48017348 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -44,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class DenseQR : QR + public class DenseQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -109,13 +108,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -154,41 +153,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs index 17845971..7c544c69 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -56,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public DenseSvd(DenseMatrix matrix, bool computeVectors) { @@ -117,13 +116,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, input.ColumnCount, dresult.Data); @@ -167,40 +166,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs new file mode 100644 index 00000000..343814ab --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Evd.cs @@ -0,0 +1,114 @@ +// +// 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.Numerics; + using Generic.Factorization; + + /// + /// Eigenvalues and eigenvectors of a real matrix. + /// + /// + /// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is + /// diagonal and the eigenvector matrix V is orthogonal. + /// I.e. A = V*D*V' and V*VT=I. + /// If A is not symmetric, then the eigenvalue matrix D is block diagonal + /// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, + /// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The + /// columns of V represent the eigenvectors in the sense that A*V = V*D, + /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly + /// conditioned, or even singular, so the validity of the equation + /// A = V*D*Inverse(V) depends upon V.Condition(). + /// + public abstract class Evd : Evd + { + /// + /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. + /// + public override double Determinant + { + get + { + var det = Complex.One; + for (var i = 0; i < VectorEv.Count; i++) + { + det *= VectorEv[i]; + + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + return 0; + } + } + + return det.Magnitude; + } + } + + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + var rank = 0; + for (var i = 0; i < VectorEv.Count; i++) + { + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + continue; + } + + rank++; + } + + return rank; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < VectorEv.Count; i++) + { + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/GramSchmidt.cs b/src/Numerics/LinearAlgebra/Double/Factorization/GramSchmidt.cs new file mode 100644 index 00000000..73eb2bac --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/GramSchmidt.cs @@ -0,0 +1,88 @@ +// +// 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 Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization. + /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix. + /// + /// + /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. + /// + public abstract class GramSchmidt : GramSchmidt + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override double Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) + { + return 0; + } + } + + return Convert.ToSingle(Math.Abs(det)); + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs new file mode 100644 index 00000000..2a99c206 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs @@ -0,0 +1,67 @@ +// +// 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 Generic.Factorization; + + /// + /// 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 : LU + { + /// + /// Gets the determinant of the matrix for which the LU factorization was computed. + /// + public override double Determinant + { + get + { + var det = 1.0; + for (var j = 0; j < Factors.RowCount; j++) + { + if (Pivots[j] != j) + { + det *= -Factors.At(j, j); + } + else + { + det *= Factors.At(j, j); + } + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs new file mode 100644 index 00000000..59fad79a --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs @@ -0,0 +1,90 @@ +// +// 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 Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition. + /// Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m) + /// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix + /// (also called right triangular matrix). + /// + /// + /// The computation of the QR decomposition is done at construction time by Householder transformation. + /// + public abstract class QR : QR + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override double Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) + { + return 0; + } + } + + return Math.Abs(det); + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/SparseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/SparseCholesky.cs deleted file mode 100644 index 4ef7ddfc..00000000 --- a/src/Numerics/LinearAlgebra/Double/Factorization/SparseCholesky.cs +++ /dev/null @@ -1,258 +0,0 @@ -// -// 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 Generic; - using Generic.Factorization; - using Properties; - - /// - /// A class which encapsulates the functionality of a Cholesky factorization for soarse matrices. - /// For a symmetric, positive definite matrix A, the Cholesky factorization - /// is an lower triangular matrix L so that A = L*L'. - /// - /// - /// 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 SparseCholesky : Cholesky - { - /// - /// Initializes a new instance of the class. This object will compute the - /// Cholesky factorization when the constructor is called and cache it's factorization. - /// - /// The matrix to factor. - /// If is null. - /// If is not a square matrix. - /// If is not positive definite. - public SparseCholesky(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). - CholeskyFactor = matrix.Clone(); - for (var j = 0; j < CholeskyFactor.RowCount; j++) - { - var d = 0.0; - for (var k = 0; k < j; k++) - { - var s = 0.0; - for (var i = 0; i < k; i++) - { - s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i); - } - - s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k); - CholeskyFactor.At(j, k, s); - d += s * s; - } - - d = matrix.At(j, j) - d; - if (d <= 0.0) - { - throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); - } - - CholeskyFactor.At(j, j, Math.Sqrt(d)); - for (var k = j + 1; k < CholeskyFactor.RowCount; k++) - { - CholeskyFactor.At(j, k, 0.0); - } - } - } - - /// - /// Solves a system of linear equations, AX = B, with A Cholesky factorized. - /// - /// The right hand side , B. - /// 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); - } - - if (result.ColumnCount != input.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); - } - - if (input.RowCount != CholeskyFactor.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - input.CopyTo(result); - var order = CholeskyFactor.RowCount; - - for (var c = 0; c < result.ColumnCount; c++) - { - // Solve L*Y = B; - double sum; - for (var i = 0; i < order; i++) - { - sum = result.At(i, c); - for (var k = i - 1; k >= 0; k--) - { - sum -= CholeskyFactor.At(i, k) * result.At(k, c); - } - - result.At(i, c, sum / CholeskyFactor.At(i, i)); - } - - // Solve L'*X = Y; - for (var i = order - 1; i >= 0; i--) - { - sum = result.At(i, c); - for (var k = i + 1; k < order; k++) - { - sum -= CholeskyFactor.At(k, i) * result.At(k, c); - } - - result.At(i, c, sum / CholeskyFactor.At(i, i)); - } - } - } - - /// - /// 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 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 != CholeskyFactor.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - input.CopyTo(result); - var order = CholeskyFactor.RowCount; - - // Solve L*Y = B; - double sum; - for (var i = 0; i < order; i++) - { - sum = result[i]; - for (var k = i - 1; k >= 0; k--) - { - sum -= CholeskyFactor.At(i, k) * result[k]; - } - - result[i] = sum / CholeskyFactor.At(i, i); - } - - // Solve L'*X = Y; - for (var i = order - 1; i >= 0; i--) - { - sum = result[i]; - for (var k = i + 1; k < order; k++) - { - sum -= CholeskyFactor.At(k, i) * result[k]; - } - - result[i] = sum / CholeskyFactor.At(i, i); - } - } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override double AddT(double val1, double val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override double LogT(double val1) - { - return Math.Log(val1); - } - #endregion - } -} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/SparseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/SparseLU.cs deleted file mode 100644 index 45d8ef9c..00000000 --- a/src/Numerics/LinearAlgebra/Double/Factorization/SparseLU.cs +++ /dev/null @@ -1,317 +0,0 @@ -// -// 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 Generic; - using Generic.Factorization; - 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 SparseLU : 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 SparseLU(Matrix 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. - var order = matrix.RowCount; - Factors = matrix.Clone(); - Pivots = new int[order]; - - // Initialize the pivot matrix to the identity permutation. - for (var i = 0; i < order; i++) - { - Pivots[i] = i; - } - - var vectorLUcolj = new double[order]; - for (var j = 0; j < order; j++) - { - // Make a copy of the j-th column to localize references. - for (var i = 0; i < order; i++) - { - vectorLUcolj[i] = Factors.At(i, j); - } - - // Apply previous transformations. - for (var i = 0; i < order; i++) - { - var kmax = Math.Min(i, j); - var s = 0.0; - for (var k = 0; k < kmax; k++) - { - s += Factors.At(i, k) * vectorLUcolj[k]; - } - - vectorLUcolj[i] -= s; - Factors.At(i, j, vectorLUcolj[i]); - } - - // Find pivot and exchange if necessary. - var p = j; - for (var i = j + 1; i < order; i++) - { - if (Math.Abs(vectorLUcolj[i]) > Math.Abs(vectorLUcolj[p])) - { - p = i; - } - } - - if (p != j) - { - for (var k = 0; k < order; k++) - { - var temp = Factors.At(p, k); - Factors.At(p, k, Factors.At(j, k)); - Factors.At(j, k, temp); - } - - Pivots[j] = p; - } - - // Compute multipliers. - if (j < order & Factors.At(j, j) != 0.0) - { - for (var i = j + 1; i < order; i++) - { - Factors.At(i, j, (Factors.At(i, j) / Factors.At(j, j))); - } - } - } - } - - /// - /// 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 != Factors.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - // Copy the contents of input to result. - input.CopyTo(result); - for (var i = 0; i < Pivots.Length; i++) - { - if (Pivots[i] == i) - { - continue; - } - - var p = Pivots[i]; - for (var j = 0; j < result.ColumnCount; j++) - { - var temp = result.At(p, j); - result.At(p, j, result.At(i, j)); - result.At(i, j, temp); - } - } - - var order = Factors.RowCount; - - // Solve L*Y = P*B - for (var k = 0; k < order; k++) - { - for (var i = k + 1; i < order; i++) - { - for (var j = 0; j < result.ColumnCount; j++) - { - var temp = result.At(k, j) * Factors.At(i, k); - result.At(i, j, result.At(i, j) - temp); - } - } - } - - // Solve U*X = Y; - for (var k = order - 1; k >= 0; k--) - { - for (var j = 0; j < result.ColumnCount; j++) - { - result.At(k, j, (result.At(k, j) / Factors.At(k, k))); - } - - for (var i = 0; i < k; i++) - { - for (var j = 0; j < result.ColumnCount; j++) - { - var temp = result.At(k, j) * Factors.At(i, k); - result.At(i, j, result.At(i, j) - temp); - } - } - } - } - - /// - /// 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 != Factors.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - // Copy the contents of input to result. - input.CopyTo(result); - for (var i = 0; i < Pivots.Length; i++) - { - if (Pivots[i] == i) - { - continue; - } - - var p = Pivots[i]; - var temp = result[p]; - result[p] = result[i]; - result[i] = temp; - } - - var order = Factors.RowCount; - - // Solve L*Y = P*B - for (var k = 0; k < order; k++) - { - for (var i = k + 1; i < order; i++) - { - result[i] -= result[k] * Factors.At(i, k); - } - } - - // Solve U*X = Y; - for (var k = order - 1; k >= 0; k--) - { - result[k] /= Factors.At(k, k); - for (var i = 0; i < k; i++) - { - result[i] -= result[k] * Factors.At(i, k); - } - } - } - - /// - /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition. - /// - /// The inverse of this matrix. - public override Matrix Inverse() - { - var order = Factors.RowCount; - var inverse = Factors.CreateMatrix(order, order); - for (var i = 0; i < order; i++) - { - inverse.At(i, i, 1.0); - } - - return Solve(inverse); - } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - #endregion - } -} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/SparseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/SparseQR.cs deleted file mode 100644 index 70a8dced..00000000 --- a/src/Numerics/LinearAlgebra/Double/Factorization/SparseQR.cs +++ /dev/null @@ -1,358 +0,0 @@ -// -// 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 System.Linq; - using Generic; - using Generic.Factorization; - using Properties; - - /// - /// A class which encapsulates the functionality of the QR decomposition. - /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal matrix - /// (its columns are orthogonal unit vectors meaning QTQ = I) and R is an upper triangular matrix - /// (also called right triangular matrix). - /// - /// - /// The computation of the QR decomposition is done at construction time by Householder transformation. - /// - public class SparseQR : QR - { - /// - /// Initializes a new instance of the class. This object will compute the - /// QR factorization when the constructor is called and cache it's factorization. - /// - /// The matrix to factor. - /// If is null. - public SparseQR(Matrix matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (matrix.RowCount < matrix.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - MatrixR = matrix.Clone(); - MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount); - - for (var i = 0; i < matrix.RowCount; i++) - { - MatrixQ.At(i, i, 1.0); - } - - var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount); - var u = new double[minmn][]; - for (var i = 0; i < minmn; i++) - { - u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i); - ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1); - } - - for (var i = minmn - 1; i >= 0; i--) - { - ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1); - } - } - - /// - /// Generate column from initial matrix to work array - /// - /// Initial matrix - /// The firts row - /// The last row - /// Column index - /// Generated vector - private static double[] GenerateColumn(Matrix a, int rowStart, int rowEnd, int column) - { - var ru = rowEnd - rowStart + 1; - var u = new double[ru]; - - for (var i = rowStart; i <= rowEnd; i++) - { - u[i - rowStart] = a.At(i, rowStart); - a.At(i, rowStart, 0.0); - } - - var norm = u.Sum(t => t * t); - norm = Math.Sqrt(norm); - - if (rowStart == rowEnd || norm == 0) - { - a.At(rowStart, column, -u[0]); - u[0] = Math.Sqrt(2.0); - return u; - } - - var scale = 1.0 / norm; - if (u[0] < 0.0) - { - scale *= -1.0; - } - - a.At(rowStart, column, -1.0 / scale); - - for (var i = 0; i < ru; i++) - { - u[i] *= scale; - } - - u[0] += 1.0; - var s = Math.Sqrt(1.0 / u[0]); - - for (var i = 0; i < ru; i++) - { - u[i] *= s; - } - - return u; - } - - /// - /// Perform calculation of Q or R - /// - /// Work array - /// Q or R matrices - /// The first row - /// The last row - /// The first column - /// The last column - private static void ComputeQR(double[] u, Matrix a, int rowStart, int rowEnd, int columnStart, int columnEnd) - { - if (rowEnd < rowStart || columnEnd < columnStart) - { - return; - } - - var v = new double[columnEnd - columnStart + 1]; - for (var j = columnStart; j <= columnEnd; j++) - { - v[j - columnStart] = 0.0; - } - - for (var i = rowStart; i <= rowEnd; i++) - { - for (var j = columnStart; j <= columnEnd; j++) - { - v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j)); - } - } - - for (var i = rowStart; i <= rowEnd; i++) - { - for (var j = columnStart; j <= columnEnd; j++) - { - a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart])); - } - } - } - - /// - /// Solves a system of linear equations, AX = B, with A QR factorized. - /// - /// The right hand side , B. - /// The left hand side , X. - public override void Solve(Matrix input, Matrix result) - { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // The solution X should have the same number of columns as B - if (input.ColumnCount != result.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); - } - - // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows - if (MatrixR.RowCount != input.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); - } - - // The solution X row dimension is equal to the column dimension of A - if (MatrixR.ColumnCount != result.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); - } - - var inputCopy = input.Clone(); - - // Compute Y = transpose(Q)*B - var bn = inputCopy.ColumnCount; - var column = new double[MatrixR.RowCount]; - for (var j = 0; j < bn; j++) - { - for (var k = 0; k < MatrixR.RowCount; k++) - { - column[k] = inputCopy.At(k, j); - } - - for (var i = 0; i < MatrixR.RowCount; i++) - { - double s = 0; - for (var k = 0; k < MatrixR.RowCount; k++) - { - s += MatrixQ.At(k, i) * column[k]; - } - - inputCopy.At(i, j, s); - } - } - - // Solve R*X = Y; - for (var k = MatrixR.ColumnCount - 1; k >= 0; k--) - { - for (var j = 0; j < bn; j++) - { - inputCopy.At(k, j, inputCopy.At(k, j) / MatrixR.At(k, k)); - } - - for (var i = 0; i < k; i++) - { - for (var j = 0; j < bn; j++) - { - inputCopy.At(i, j, inputCopy.At(i, j) - (inputCopy.At(k, j) * MatrixR.At(i, k))); - } - } - } - - for (var i = 0; i < MatrixR.ColumnCount; i++) - { - for (var j = 0; j < inputCopy.ColumnCount; j++) - { - result.At(i, j, inputCopy.At(i, j)); - } - } - } - - /// - /// Solves a system of linear equations, Ax = b, with A QR factorized. - /// - /// The right hand side vector, b. - /// The left hand side , x. - public override void Solve(Vector input, Vector result) - { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Ax=b where A is an m x n matrix - // Check that b is a column vector with m entries - if (MatrixR.RowCount != input.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength); - } - - // Check that x is a column vector with n entries - if (MatrixR.ColumnCount != result.Count) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - var inputCopy = input.Clone(); - - // Compute Y = transpose(Q)*B - var column = new double[MatrixR.RowCount]; - for (var k = 0; k < MatrixR.RowCount; k++) - { - column[k] = inputCopy[k]; - } - - for (var i = 0; i < MatrixR.RowCount; i++) - { - double s = 0; - for (var k = 0; k < MatrixR.RowCount; k++) - { - s += MatrixQ.At(k, i) * column[k]; - } - - inputCopy[i] = s; - } - - // Solve R*X = Y; - for (var k = MatrixR.ColumnCount - 1; k >= 0; k--) - { - inputCopy[k] /= MatrixR.At(k, k); - for (var i = 0; i < k; i++) - { - inputCopy[i] -= inputCopy[k] * MatrixR.At(i, k); - } - } - - for (var i = 0; i < MatrixR.ColumnCount; i++) - { - result[i] = inputCopy[i]; - } - } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - #endregion - } -} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/SparseSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/SparseSvd.cs deleted file mode 100644 index 00c87a56..00000000 --- a/src/Numerics/LinearAlgebra/Double/Factorization/SparseSvd.cs +++ /dev/null @@ -1,950 +0,0 @@ -// -// 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 Generic; - using Generic.Factorization; - using Properties; - - /// - /// A class which encapsulates the functionality of the singular value decomposition (SVD) for . - /// Suppose M is an m-by-n matrix whose entries are real numbers. - /// Then there exists a factorization of the form M = UΣVT where: - /// - U is an m-by-m unitary matrix; - /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; - /// - VT denotes transpose of V, an n-by-n unitary matrix; - /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal - /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined - /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M. - /// - /// - /// The computation of the singular value decomposition is done at construction time. - /// - public class SparseSvd : Svd - { - /// - /// Initializes a new instance of the class. This object will compute the - /// the singular value decomposition when the constructor is called and cache it's decomposition. - /// - /// The matrix to factor. - /// Compute the singular U and VT vectors or not. - /// If is null. - /// If SVD algorithm failed to converge with matrix . - public SparseSvd(Matrix matrix, bool computeVectors) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - ComputeVectors = computeVectors; - var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); - var matrixCopy = matrix.Clone(); - - VectorS = matrixCopy.CreateVector(nm); - MatrixU = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); - MatrixVT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); - - const int Maxiter = 1000; - var e = new double[matrixCopy.ColumnCount]; - var work = new double[matrixCopy.RowCount]; - - int i, j; - int l, lp1; - var cs = 0.0; - var sn = 0.0; - double t; - - var ncu = matrixCopy.RowCount; - - // Reduce matrixCopy to bidiagonal form, storing the diagonal elements - // In s and the super-diagonal elements in e. - var nct = Math.Min(matrixCopy.RowCount - 1, matrixCopy.ColumnCount); - var nrt = Math.Max(0, Math.Min(matrixCopy.ColumnCount - 2, matrixCopy.RowCount)); - var lu = Math.Max(nct, nrt); - for (l = 0; l < lu; l++) - { - lp1 = l + 1; - if (l < nct) - { - // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l]. - var xnorm = Dnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); - VectorS[l] = xnorm; - if (VectorS[l] != 0.0) - { - if (matrixCopy.At(l, l) != 0.0) - { - VectorS[l] = Dsign(VectorS[l], matrixCopy.At(l, l)); - } - - DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0 / VectorS[l]); - matrixCopy.At(l, l, (1.0 + matrixCopy.At(l, l))); - } - - VectorS[l] = -VectorS[l]; - } - - for (j = lp1; j < matrixCopy.ColumnCount; j++) - { - if (l < nct) - { - if (VectorS[l] != 0.0) - { - // Apply the transformation. - t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); - for (var ii = l; ii < matrixCopy.RowCount; ii++) - { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); - } - } - } - - // Place the l-th row of matrixCopy into e for the - // Subsequent calculation of the row transformation. - e[j] = matrixCopy.At(l, j); - } - - if (ComputeVectors && l < nct) - { - // Place the transformation in u for subsequent back multiplication. - for (i = l; i < matrixCopy.RowCount; i++) - { - MatrixU.At(i, l, matrixCopy.At(i, l)); - } - } - - if (l >= nrt) - { - continue; - } - - // Compute the l-th row transformation and place the l-th super-diagonal in e(l). - var enorm = Dnrm2Vector(e, lp1); - e[l] = enorm; - if (e[l] != 0.0) - { - if (e[lp1] != 0.0) - { - e[l] = Dsign(e[l], e[lp1]); - } - - DscalVector(e, lp1, 1.0 / e[l]); - e[lp1] = 1.0 + e[lp1]; - } - - e[l] = -e[l]; - if (lp1 < matrixCopy.RowCount && e[l] != 0.0) - { - // Apply the transformation. - for (i = lp1; i < matrixCopy.RowCount; i++) - { - work[i] = 0.0; - } - - for (j = lp1; j < matrixCopy.ColumnCount; j++) - { - for (var ii = lp1; ii < matrixCopy.RowCount; ii++) - { - work[ii] += e[j] * matrixCopy.At(ii, j); - } - } - - for (j = lp1; j < matrixCopy.ColumnCount; j++) - { - var ww = -e[j] / e[lp1]; - for (var ii = lp1; ii < matrixCopy.RowCount; ii++) - { - matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); - } - } - } - - if (ComputeVectors) - { - // Place the transformation in v for subsequent back multiplication. - for (i = lp1; i < matrixCopy.ColumnCount; i++) - { - MatrixVT.At(i, l, e[i]); - } - } - } - - // Set up the final bidiagonal matrixCopy or order m. - var m = Math.Min(matrixCopy.ColumnCount, matrixCopy.RowCount + 1); - var nctp1 = nct + 1; - var nrtp1 = nrt + 1; - if (nct < matrixCopy.ColumnCount) - { - VectorS[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); - } - - if (matrixCopy.RowCount < m) - { - VectorS[m - 1] = 0.0; - } - - if (nrtp1 < m) - { - e[nrtp1 - 1] = matrixCopy.At((nrtp1 - 1), (m - 1)); - } - - e[m - 1] = 0.0; - - // If required, generate u. - if (ComputeVectors) - { - for (j = nctp1 - 1; j < ncu; j++) - { - for (i = 0; i < matrixCopy.RowCount; i++) - { - MatrixU.At(i, j, 0.0); - } - - MatrixU.At(j, j, 1.0); - } - - for (l = nct - 1; l >= 0; l--) - { - if (VectorS[l] != 0.0) - { - for (j = l + 1; j < ncu; j++) - { - t = -Ddot(MatrixU, matrixCopy.RowCount, l, j, l) / MatrixU.At(l, l); - for (var ii = l; ii < matrixCopy.RowCount; ii++) - { - MatrixU.At(ii, j, MatrixU.At(ii, j) + (t * MatrixU.At(ii, l))); - } - } - - DscalColumn(MatrixU, matrixCopy.RowCount, l, l, -1.0); - MatrixU.At(l, l, 1.0 + MatrixU.At(l, l)); - for (i = 0; i < l; i++) - { - MatrixU.At(i, l, 0.0); - } - } - else - { - for (i = 0; i < matrixCopy.RowCount; i++) - { - MatrixU.At(i, l, 0.0); - } - - MatrixU.At(l, l, 1.0); - } - } - } - - // If it is required, generate v. - if (ComputeVectors) - { - for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) - { - lp1 = l + 1; - if (l < nrt) - { - if (e[l] != 0.0) - { - for (j = lp1; j < matrixCopy.ColumnCount; j++) - { - t = -Ddot(MatrixVT, matrixCopy.ColumnCount, l, j, lp1) / MatrixVT.At(lp1, l); - for (var ii = l; ii < matrixCopy.ColumnCount; ii++) - { - MatrixVT.At(ii, j, MatrixVT.At(ii, j) + (t * MatrixVT.At(ii, l))); - } - } - } - } - - for (i = 0; i < matrixCopy.ColumnCount; i++) - { - MatrixVT.At(i, l, 0.0); - } - - MatrixVT.At(l, l, 1.0); - } - } - - // Transform s and e so that they are double . - for (i = 0; i < m; i++) - { - double r; - if (VectorS[i] != 0.0) - { - t = VectorS[i]; - r = VectorS[i] / t; - VectorS[i] = t; - if (i < m - 1) - { - e[i] = e[i] / r; - } - - if (ComputeVectors) - { - DscalColumn(MatrixU, matrixCopy.RowCount, i, 0, r); - } - } - - // Exit - if (i == m - 1) - { - break; - } - - if (e[i] != 0.0) - { - t = e[i]; - r = t / e[i]; - e[i] = t; - VectorS[i + 1] = VectorS[i + 1] * r; - if (ComputeVectors) - { - DscalColumn(MatrixVT, matrixCopy.ColumnCount, i + 1, 0, r); - } - } - } - - // Main iteration loop for the singular values. - var mn = m; - var iter = 0; - - while (m > 0) - { - // Quit if all the singular values have been found. If too many iterations have been performed, - // throw exception that Convergence Failed - if (iter >= Maxiter) - { - throw new ArgumentException(Resources.ConvergenceFailed); - } - - // This section of the program inspects for negligible elements in the s and e arrays. On - // completion the variables kase and l are set as follows. - // Kase = 1 if VectorS[m] and e[l-1] are negligible and l < m - // Kase = 2 if VectorS[l] is negligible and l < m - // Kase = 3 if e[l-1] is negligible, l < m, and VectorS[l, ..., VectorS[m] are not negligible (qr step). - // Лase = 4 if e[m-1] is negligible (convergence). - double ztest; - double test; - for (l = m - 2; l >= 0; l--) - { - test = Math.Abs(VectorS[l]) + Math.Abs(VectorS[l + 1]); - ztest = test + Math.Abs(e[l]); - if (ztest.AlmostEqualInDecimalPlaces(test, 15)) - { - e[l] = 0.0; - break; - } - } - - int kase; - if (l == m - 2) - { - kase = 4; - } - else - { - int ls; - for (ls = m - 1; ls > l; ls--) - { - test = 0.0; - if (ls != m - 1) - { - test = test + Math.Abs(e[ls]); - } - - if (ls != l + 1) - { - test = test + Math.Abs(e[ls - 1]); - } - - ztest = test + Math.Abs(VectorS[ls]); - if (ztest.AlmostEqualInDecimalPlaces(test, 15)) - { - VectorS[ls] = 0.0; - break; - } - } - - if (ls == l) - { - kase = 3; - } - else if (ls == m - 1) - { - kase = 1; - } - else - { - kase = 2; - l = ls; - } - } - - l = l + 1; - - // Perform the task indicated by kase. - int k; - double f; - switch (kase) - { - // Deflate negligible VectorS[m]. - case 1: - f = e[m - 2]; - e[m - 2] = 0.0; - double t1; - for (var kk = l; kk < m - 1; kk++) - { - k = m - 2 - kk + l; - t1 = VectorS[k]; - Drotg(ref t1, ref f, ref cs, ref sn); - VectorS[k] = t1; - if (k != l) - { - f = -sn * e[k - 1]; - e[k - 1] = cs * e[k - 1]; - } - - if (ComputeVectors) - { - Drot(MatrixVT, matrixCopy.ColumnCount, k, m - 1, cs, sn); - } - } - - break; - - // Split at negligible VectorS[l]. - case 2: - f = e[l - 1]; - e[l - 1] = 0.0; - for (k = l; k < m; k++) - { - t1 = VectorS[k]; - Drotg(ref t1, ref f, ref cs, ref sn); - VectorS[k] = t1; - f = -sn * e[k]; - e[k] = cs * e[k]; - if (ComputeVectors) - { - Drot(MatrixU, matrixCopy.RowCount, k, l - 1, cs, sn); - } - } - - break; - - // Perform one qr step. - case 3: - // Calculate the shift. - var scale = 0.0; - scale = Math.Max(scale, Math.Abs(VectorS[m - 1])); - scale = Math.Max(scale, Math.Abs(VectorS[m - 2])); - scale = Math.Max(scale, Math.Abs(e[m - 2])); - scale = Math.Max(scale, Math.Abs(VectorS[l])); - scale = Math.Max(scale, Math.Abs(e[l])); - var sm = VectorS[m - 1] / scale; - var smm1 = VectorS[m - 2] / scale; - var emm1 = e[m - 2] / scale; - var sl = VectorS[l] / scale; - var el = e[l] / scale; - var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0; - var c = (sm * emm1) * (sm * emm1); - var shift = 0.0; - if (b != 0.0 || c != 0.0) - { - shift = Math.Sqrt((b * b) + c); - if (b < 0.0) - { - shift = -shift; - } - - shift = c / (b + shift); - } - - f = ((sl + sm) * (sl - sm)) + shift; - var g = sl * el; - - // Chase zeros. - for (k = l; k < m - 1; k++) - { - Drotg(ref f, ref g, ref cs, ref sn); - if (k != l) - { - e[k - 1] = f; - } - - f = (cs * VectorS[k]) + (sn * e[k]); - e[k] = (cs * e[k]) - (sn * VectorS[k]); - g = sn * VectorS[k + 1]; - VectorS[k + 1] = cs * VectorS[k + 1]; - if (ComputeVectors) - { - Drot(MatrixVT, matrixCopy.ColumnCount, k, k + 1, cs, sn); - } - - Drotg(ref f, ref g, ref cs, ref sn); - VectorS[k] = f; - f = (cs * e[k]) + (sn * VectorS[k + 1]); - VectorS[k + 1] = (-sn * e[k]) + (cs * VectorS[k + 1]); - g = sn * e[k + 1]; - e[k + 1] = cs * e[k + 1]; - if (ComputeVectors && k < matrixCopy.RowCount) - { - Drot(MatrixU, matrixCopy.RowCount, k, k + 1, cs, sn); - } - } - - e[m - 2] = f; - iter = iter + 1; - break; - - // Convergence. - case 4: - // Make the singular value positive - if (VectorS[l] < 0.0) - { - VectorS[l] = -VectorS[l]; - if (ComputeVectors) - { - DscalColumn(MatrixVT, matrixCopy.ColumnCount, l, 0, -1.0); - } - } - - // Order the singular value. - while (l != mn - 1) - { - if (VectorS[l] >= VectorS[l + 1]) - { - break; - } - - t = VectorS[l]; - VectorS[l] = VectorS[l + 1]; - VectorS[l + 1] = t; - if (ComputeVectors && l < matrixCopy.ColumnCount) - { - Dswap(MatrixVT, matrixCopy.ColumnCount, l, l + 1); - } - - if (ComputeVectors && l < matrixCopy.RowCount) - { - Dswap(MatrixU, matrixCopy.RowCount, l, l + 1); - } - - l = l + 1; - } - - iter = 0; - m = m - 1; - break; - } - } - - if (ComputeVectors) - { - MatrixVT = MatrixVT.Transpose(); - } - - // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses - // a singular vector of length mRows+1 when mRows < mColumns. The last element is not used and needs to be removed. - // we should port lapack's svd routine to remove this problem. - if (matrixCopy.RowCount < matrixCopy.ColumnCount) - { - nm--; - var tmp = matrixCopy.CreateVector(nm); - for (i = 0; i < nm; i++) - { - tmp[i] = VectorS[i]; - } - - VectorS = tmp; - } - } - - /// - /// Calculates absolute value of multiplied on signum function of - /// - /// Double value z1 - /// Double value z2 - /// Result multiplication of signum function and absolute value - private static double Dsign(double z1, double z2) - { - return Math.Abs(z1) * (z2 / Math.Abs(z2)); - } - - /// - /// Swap column and - /// - /// Source matrix - /// The number of rows in - /// Column A index to swap - /// Column B index to swap - private static void Dswap(Matrix a, int rowCount, int columnA, int columnB) - { - for (var i = 0; i < rowCount; i++) - { - var z = a.At(i, columnA); - a.At(i, columnA, a.At(i, columnB)); - a.At(i, columnB, z); - } - } - - /// - /// Scale column by starting from row - /// - /// Source matrix - /// The number of rows in - /// Column to scale - /// Row to scale from - /// Scale value - private static void DscalColumn(Matrix a, int rowCount, int column, int rowStart, double z) - { - for (var i = rowStart; i < rowCount; i++) - { - a.At(i, column, a.At(i, column) * z); - } - } - - /// - /// Scale vector by starting from index - /// - /// Source vector - /// Row to scale from - /// Scale value - private static void DscalVector(double[] a, int start, double z) - { - for (var i = start; i < a.Length; i++) - { - a[i] = a[i] * z; - } - } - - /// - /// Given the Cartesian coordinates (da, db) of a point p, these fucntion 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 - /// Provides the y-coordinate of the point p. On exit contains the parameter z associated with the Givens rotation - /// Contains the parameter c associated with the Givens rotation - /// Contains the parameter s associated with the Givens rotation - /// This is equivalent to the DROTG LAPACK routine. - private static void Drotg(ref double da, ref double db, ref double c, ref double s) - { - double r, z; - - var roe = db; - var absda = Math.Abs(da); - var absdb = Math.Abs(db); - if (absda > absdb) - { - roe = da; - } - - var scale = absda + absdb; - if (scale == 0.0) - { - c = 1.0; - s = 0.0; - r = 0.0; - z = 0.0; - } - else - { - var sda = da / scale; - var sdb = db / scale; - r = scale * Math.Sqrt((sda * sda) + (sdb * sdb)); - if (roe < 0.0) - { - r = -r; - } - - c = da / r; - s = db / r; - z = 1.0; - if (absda > absdb) - { - z = s; - } - - if (absdb >= absda && c != 0.0) - { - z = 1.0 / c; - } - } - - da = r; - db = z; - } - - /// dded - /// Calculate Norm 2 of the column in matrix starting from row - /// - /// Source matrix - /// The number of rows in - /// Column index - /// Start row index - /// Norm2 (Euclidean norm) of trhe column - private static double Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) - { - double s = 0; - for (var i = rowStart; i < rowCount; i++) - { - s += a.At(i, column) * a.At(i, column); - } - - return Math.Sqrt(s); - } - - /// - /// Calculate Norm 2 of the vector starting from index - /// - /// Source vector - /// Start index - /// Norm2 (Euclidean norm) of the vector - private static double Dnrm2Vector(double[] a, int rowStart) - { - double s = 0; - for (var i = rowStart; i < a.Length; i++) - { - s += a[i] * a[i]; - } - - return Math.Sqrt(s); - } - - /// - /// Calculate dot product of and - /// - /// Source matrix - /// The number of rows in - /// Index of column A - /// Index of column B - /// Starting row index - /// Dot product value - private static double Ddot(Matrix a, int rowCount, int columnA, int columnB, int rowStart) - { - var z = 0.0; - for (var i = rowStart; i < rowCount; i++) - { - z += a.At(i, columnB) * a.At(i, columnA); - } - - return z; - } - - /// - /// Performs rotation of points in the plane. Given two vectors x and y , - /// each vector element of these vectors is replaced as follows: x(i) = c*x(i) + s*y(i); y(i) = c*y(i) - s*x(i) - /// - /// Source matrix - /// The number of rows in - /// Index of column A - /// Index of column B - /// Scalar "c" value - /// Scalar "s" value - private static void Drot(Matrix a, int rowCount, int columnA, int columnB, double c, double s) - { - for (var i = 0; i < rowCount; i++) - { - var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); - var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); - a.At(i, columnB, tmp); - a.At(i, columnA, z); - } - } - - /// - /// Solves a system of linear equations, AX = B, with A SVD 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"); - } - - if (!ComputeVectors) - { - throw new InvalidOperationException(Resources.SingularVectorsNotComputed); - } - - // The solution X should have the same number of columns as B - if (input.ColumnCount != result.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); - } - - // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows - if (MatrixU.RowCount != input.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); - } - - // The solution X row dimension is equal to the column dimension of A - if (MatrixVT.ColumnCount != result.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); - } - - var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); - var bn = input.ColumnCount; - - var tmp = new double[MatrixVT.ColumnCount]; - - for (var k = 0; k < bn; k++) - { - for (var j = 0; j < MatrixVT.ColumnCount; j++) - { - double value = 0; - if (j < mn) - { - for (var i = 0; i < MatrixU.RowCount; i++) - { - value += MatrixU.At(i, j) * input.At(i, k); - } - - value /= VectorS[j]; - } - - tmp[j] = value; - } - - for (var j = 0; j < MatrixVT.ColumnCount; j++) - { - double value = 0; - for (var i = 0; i < MatrixVT.ColumnCount; i++) - { - value += MatrixVT.At(i, j) * tmp[i]; - } - - result[j, k] = value; - } - } - } - - /// - /// Solves a system of linear equations, Ax = b, with A SVD factorized. - /// - /// The right hand side vector, b. - /// The left hand side , x. - public override void Solve(Vector input, Vector result) - { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (!ComputeVectors) - { - throw new InvalidOperationException(Resources.SingularVectorsNotComputed); - } - - // Ax=b where A is an m x n matrix - // Check that b is a column vector with m entries - if (MatrixU.RowCount != input.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength); - } - - // Check that x is a column vector with n entries - if (MatrixVT.ColumnCount != result.Count) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); - var tmp = new double[MatrixVT.ColumnCount]; - double value; - for (var j = 0; j < MatrixVT.ColumnCount; j++) - { - value = 0; - if (j < mn) - { - for (var i = 0; i < MatrixU.RowCount; i++) - { - value += MatrixU.At(i, j) * input[i]; - } - - value /= VectorS[j]; - } - - tmp[j] = value; - } - - for (var j = 0; j < MatrixVT.ColumnCount; j++) - { - value = 0; - for (int i = 0; i < MatrixVT.ColumnCount; i++) - { - value += MatrixVT.At(i, j) * tmp[i]; - } - - result[j] = value; - } - } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - - #endregion - } -} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs new file mode 100644 index 00000000..238bda49 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs @@ -0,0 +1,118 @@ +// +// 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 System.Linq; + using Generic; + using Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the singular value decomposition (SVD). + /// Suppose M is an m-by-n matrix whose entries are real numbers. + /// Then there exists a factorization of the form M = UΣVT where: + /// - U is an m-by-m unitary matrix; + /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; + /// - VT denotes transpose of V, an n-by-n unitary matrix; + /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal + /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined + /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M. + /// + /// + /// The computation of the singular value decomposition is done at construction time. + /// + public abstract class Svd : Svd + { + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + return VectorS.Count(t => !Math.Abs(t).AlmostEqual(0.0)); + } + } + + /// + /// Gets the two norm of the . + /// + /// The 2-norm of the . + public override double Norm2 + { + get + { + return Math.Abs(VectorS[0]); + } + } + + /// + /// Gets the condition number max(S) / min(S) + /// + /// The condition number. + public override double ConditionNumber + { + get + { + var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; + return Math.Abs(VectorS[0]) / Math.Abs(VectorS[tmp]); + } + } + + /// + /// Gets the determinant of the square matrix for which the SVD was computed. + /// + public override double Determinant + { + get + { + if (MatrixU.RowCount != MatrixVT.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + foreach (var value in VectorS) + { + det *= value; + if (Math.Abs(value).AlmostEqual(0.0)) + { + return 0; + } + } + + return Math.Abs(det); + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs index 91dacd86..28c719f1 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -44,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 class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -220,39 +219,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization result[i] = sum / CholeskyFactor.At(i, i); } } - - #region Simple T Mathematics - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override double AddT(double val1, double val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override double LogT(double val1) - { - return Math.Log(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs index 3850d09f..a9027250 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,16 +47,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public UserEvd(Matrix matrix) { @@ -1218,16 +1217,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserGramSchmidt.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserGramSchmidt.cs index 417a5922..81b91085 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserGramSchmidt.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -42,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class UserGramSchmidt : GramSchmidt + public class UserGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an orthogonal matrix @@ -244,29 +243,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs index e29c68cb..d6e1701b 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class UserLU : LU + public class UserLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -298,20 +297,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization return Solve(inverse); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs index f3f349c2..c9d6695c 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization using System; using System.Linq; using Generic; - using Generic.Factorization; using Properties; /// @@ -45,7 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class UserQR : QR + public class UserQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -91,7 +90,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Generate column from initial matrix to work array /// /// Initial matrix - /// The firts row + /// The first row /// The last row /// Column index /// Generated vector @@ -329,29 +328,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs index 1eed31dd..96f82f33 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -56,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public UserSvd(Matrix matrix, bool computeVectors) { @@ -702,14 +701,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization db = z; } - /// dded + /// /// Calculate Norm 2 of the column in matrix starting from row /// /// Source matrix /// The number of rows in /// Column index /// Start row index - /// Norm2 (Euclidean norm) of trhe column + /// Norm2 (Euclidean norm) of the column private static double Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) { double s = 0; @@ -921,30 +920,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization result[j] = value; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/Cholesky.cs index fa205220..562cc87b 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/Cholesky.cs @@ -99,7 +99,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization return new LinearAlgebra.Complex32.Factorization.UserCholesky(matrix as Matrix) as Cholesky; } - throw new NotImplementedException(); + throw new NotSupportedException(); } /// @@ -125,36 +125,17 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// /// Gets the determinant of the matrix for which the Cholesky matrix was computed. /// - public virtual T Determinant + public abstract T Determinant { - get - { - var det = OneValueT; - for (var j = 0; j < CholeskyFactor.RowCount; j++) - { - det = MultiplyT(det, MultiplyT(CholeskyFactor[j, j], CholeskyFactor[j, j])); - } - - return det; - } + get; } /// /// Gets the log determinant of the matrix for which the Cholesky matrix was computed. /// - public virtual T DeterminantLn + public abstract T DeterminantLn { - get - { - var det = default(T); - for (var j = 0; j < CholeskyFactor.RowCount; j++) - { - // det += 2.0 * CholeskyFactor[j, j].NaturalLogarithm(); - det = AddT(det, MultiplyT(AddT(OneValueT, OneValueT), LogT(CholeskyFactor[j, j]))); - } - - return det; - } + get; } /// @@ -206,67 +187,5 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// The right hand side vector, b. /// The left hand side , x. public abstract void Solve(Vector input, Vector result); - - #region Simple arithmetic of type T - - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected abstract T AddT(T val1, T val2); - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected abstract T MultiplyT(T val1, T val2); - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected abstract T LogT(T val1); - - /// - /// Gets value of type T equal to one - /// - /// One value - private static T OneValueT - { - get - { - if (typeof(T) == typeof(Complex)) - { - object one = Complex.One; - return (T)one; - } - - if (typeof(T) == typeof(Complex32)) - { - object one = Complex32.One; - return (T)one; - } - - if (typeof(T) == typeof(double)) - { - object one = 1.0d; - return (T)one; - } - - if (typeof(T) == typeof(float)) - { - object one = 1.0f; - return (T)one; - } - - throw new NotSupportedException(); - } - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/Evd.cs index faff97fb..b14622b2 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/Evd.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/Evd.cs @@ -48,7 +48,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// /// Supported data types are double, single, , and . public abstract class Evd : ISolver @@ -63,6 +63,32 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization protected set; } + /// + /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. + /// + public abstract T Determinant + { + get; + } + + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public abstract int Rank + { + get; + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public abstract bool IsFullRank + { + get; + } + /// /// Gets or sets the eigen values (λ) of matrix in ascending value. /// @@ -141,104 +167,19 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization return new LinearAlgebra.Complex32.Factorization.UserEvd(matrix as Matrix) as Evd; } - throw new NotImplementedException(); - } - - /// - /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. - /// - public virtual double Determinant - { - get - { - var det = Complex.One; - for (var i = 0; i < VectorEv.Count; i++) - { - det *= VectorEv[i]; - - if (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) - { - if (((Complex32)VectorEv[i]).AlmostEqual(Complex32.Zero)) - { - return 0; - } - } - else - { - if (VectorEv[i].AlmostEqual(Complex.Zero)) - { - return 0; - } - } - } - - return det.Magnitude; - } - } - - /// - /// Gets the effective numerical matrix rank. - /// - /// The number of non-negligible singular values. - public virtual int Rank - { - get - { - var rank = 0; - for (var i = 0; i < VectorEv.Count; i++) - { - if (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) - { - if (((Complex32)VectorEv[i]).AlmostEqual(Complex32.Zero)) - { - continue; - } - } - else - { - if (VectorEv[i].AlmostEqual(Complex.Zero)) - { - continue; - } - } - - rank++; - } - - return rank; - } - } - - /// - /// Gets a value indicating whether the matrix is full rank or not. - /// - /// true if the matrix is full rank; otherwise false. - public virtual bool IsFullRank - { - get - { - for (var i = 0; i < VectorEv.Count; i++) - { - if (VectorEv[i].AlmostEqual(Complex.Zero)) - { - return false; - } - } - - return true; - } + throw new NotSupportedException(); } /// Returns the eigen values as a . /// The eigen values. - public Vector EValues() + public Vector EigenValues() { return VectorEv.Clone(); } /// Returns the right eigen vectors as a . /// The eigen vectors. - public Matrix EVectors() + public Matrix EigenVectors() { return MatrixEv.Clone(); } @@ -299,16 +240,5 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// The right hand side vector, b. /// The left hand side , x. public abstract void Solve(Vector input, Vector result); - - #region Simple arithmetic of type T - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected abstract T MultiplyT(T val1, T val2); - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/GramSchmidt.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/GramSchmidt.cs index 94b43bd7..70bb4722 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/GramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/GramSchmidt.cs @@ -30,7 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization using System.Numerics; using Generic; using Numerics; - using Properties; /// /// A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization. @@ -94,45 +93,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization return new LinearAlgebra.Complex32.Factorization.UserGramSchmidt(matrix as Matrix) as GramSchmidt; } - throw new NotImplementedException(); - } - - /// - /// Gets a value indicating whether the matrix is full rank or not. - /// - /// true if the matrix is full rank; otherwise false. - public sealed override bool IsFullRank - { - get - { - return true; - } - } - - /// - /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. - /// - public override double Determinant - { - get - { - if (MatrixQ.RowCount != MatrixQ.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSquare); - } - - var det = OneValueT; - for (var i = 0; i < MatrixR.ColumnCount; i++) - { - det = MultiplyT(det, MatrixR.At(i, i)); - if (AbsoluteT(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) ? 7 : 15)) - { - return 0; - } - } - - return AbsoluteT(det); - } + throw new NotSupportedException(); } } } diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/LU.cs index 11416ea5..2118ade4 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/LU.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/LU.cs @@ -45,6 +45,11 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization public abstract class LU : ISolver where T : struct, IEquatable, IFormattable { + /// + /// Value of one for T. + /// + private static readonly T One = Common.SetOne(); + /// /// Gets or sets both the L and U factors in the same matrix. /// @@ -114,7 +119,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization return new LinearAlgebra.Complex32.Factorization.UserLU(matrix as Matrix) as LU; } - throw new NotImplementedException(); + throw new NotSupportedException(); } /// @@ -127,7 +132,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization var result = Factors.LowerTriangle(); for (var i = 0; i < result.RowCount; i++) { - result.At(i, i, OneValueT); + result.At(i, i, One); } return result; @@ -159,25 +164,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// /// Gets the determinant of the matrix for which the LU factorization was computed. /// - public virtual T Determinant + public abstract T Determinant { - get - { - var det = OneValueT; - for (var j = 0; j < Factors.RowCount; j++) - { - if (Pivots[j] != j) - { - det = MultiplyT(MinusOneValueT, MultiplyT(det, Factors.At(j, j))); - } - else - { - det = MultiplyT(det, Factors.At(j, j)); - } - } - - return det; - } + get; } /// @@ -235,89 +224,5 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// /// The inverse of this matrix. public abstract Matrix Inverse(); - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected abstract T MultiplyT(T val1, T val2); - - /// - /// Gets value of type T equal to one - /// - /// One value - private static T OneValueT - { - get - { - if (typeof(T) == typeof(Complex)) - { - object one = Complex.One; - return (T)one; - } - - if (typeof(T) == typeof(Complex32)) - { - object one = Complex32.One; - return (T)one; - } - - if (typeof(T) == typeof(double)) - { - object one = 1.0d; - return (T)one; - } - - if (typeof(T) == typeof(float)) - { - object one = 1.0f; - return (T)one; - } - - throw new NotSupportedException(); - } - } - - /// - /// Gets value of type T equal to one - /// - /// One value - private static T MinusOneValueT - { - get - { - if (typeof(T) == typeof(Complex)) - { - object one = -Complex.One; - return (T)one; - } - - if (typeof(T) == typeof(Complex32)) - { - object one = -Complex32.One; - return (T)one; - } - - if (typeof(T) == typeof(double)) - { - object one = -1.0d; - return (T)one; - } - - if (typeof(T) == typeof(float)) - { - object one = -1.0f; - return (T)one; - } - - throw new NotSupportedException(); - } - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs index cc36babe..ea4c65e3 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs @@ -114,7 +114,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization return new LinearAlgebra.Complex32.Factorization.UserQR(matrix as Matrix) as QR; } - throw new NotImplementedException(); + throw new NotSupportedException(); } /// @@ -142,47 +142,18 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. /// - public virtual double Determinant + public abstract T Determinant { - get - { - if (MatrixR.RowCount != MatrixR.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSquare); - } - - var det = OneValueT; - for (var i = 0; i < MatrixR.ColumnCount; i++) - { - det = MultiplyT(det, MatrixR.At(i, i)); - if (AbsoluteT(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) ? 7 : 15)) - { - return 0; - } - } - - return AbsoluteT(det); - } + get; } /// /// Gets a value indicating whether the matrix is full rank or not. /// /// true if the matrix is full rank; otherwise false. - public virtual bool IsFullRank + public abstract bool IsFullRank { - get - { - for (var i = 0; i < MatrixR.ColumnCount; i++) - { - if (AbsoluteT(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) ? 7 : 15)) - { - return false; - } - } - - return true; - } + get; } /// @@ -234,59 +205,5 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// The right hand side vector, b. /// The left hand side , x. public abstract void Solve(Vector input, Vector result); - - #region Simple arithmetic of type T - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected abstract T MultiplyT(T val1, T val2); - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected abstract double AbsoluteT(T val); - - /// - /// Gets value of type T equal to one - /// - /// One value - protected static T OneValueT - { - get - { - if (typeof(T) == typeof(Complex)) - { - object one = Complex.One; - return (T)one; - } - - if (typeof(T) == typeof(Complex32)) - { - object one = Complex32.One; - return (T)one; - } - - if (typeof(T) == typeof(double)) - { - object one = 1.0d; - return (T)one; - } - - if (typeof(T) == typeof(float)) - { - object one = 1.0f; - return (T)one; - } - - throw new NotSupportedException(); - } - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs index 1ba48da3..fae74399 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs @@ -95,12 +95,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// Gets the effective numerical matrix rank. /// /// The number of non-negligible singular values. - public virtual int Rank + public abstract int Rank { - get - { - return VectorS.Count(t => !AbsoluteT(t).AlmostEqualInDecimalPlaces(0.0, (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) ? 7 : 15)); - } + get; } /// @@ -155,59 +152,33 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization return new LinearAlgebra.Complex32.Factorization.UserSvd(matrix as Matrix, computeVectors) as Svd; } - throw new NotImplementedException(); + throw new NotSupportedException(); } /// /// Gets the two norm of the . /// /// The 2-norm of the . - public virtual T Norm2 + public abstract T Norm2 { - get - { - throw new NotImplementedException(); - //return AbsoluteT(VectorS[0]); - } + get; } /// /// Gets the condition number max(S) / min(S) /// /// The condition number. - public virtual double ConditionNumber + public abstract T ConditionNumber { - get - { - var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; - return AbsoluteT(VectorS[0]) / AbsoluteT(VectorS[tmp]); - } + get; } /// /// Gets the determinant of the square matrix for which the SVD was computed. /// - public virtual double Determinant + public abstract T Determinant { - get - { - if (MatrixU.RowCount != MatrixVT.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSquare); - } - - var det = OneValueT; - for (var i = 0; i < VectorS.Count; i++) - { - det = MultiplyT(det, VectorS[i]); - if (AbsoluteT(VectorS[i]).AlmostEqualInDecimalPlaces(0.0, (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) ? 7 : 15)) - { - return 0; - } - } - - return AbsoluteT(det); - } + get; } /// Returns the left singular vectors as a . @@ -312,59 +283,5 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// The right hand side vector, b. /// The left hand side , x. public abstract void Solve(Vector input, Vector result); - - #region Simple arithmetic of type T - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected abstract T MultiplyT(T val1, T val2); - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected abstract double AbsoluteT(T val); - - /// - /// Gets value of type T equal to one - /// - /// One value - private static T OneValueT - { - get - { - if (typeof(T) == typeof(Complex)) - { - object one = Complex.One; - return (T)one; - } - - if (typeof(T) == typeof(Complex32)) - { - object one = Complex32.One; - return (T)one; - } - - if (typeof(T) == typeof(double)) - { - object one = 1.0d; - return (T)one; - } - - if (typeof(T) == typeof(float)) - { - object one = 1.0f; - return (T)one; - } - - throw new NotSupportedException(); - } - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs new file mode 100644 index 00000000..6eddd377 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs @@ -0,0 +1,81 @@ +// +// 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.Single.Factorization +{ + using System; + using Generic.Factorization; + + /// + /// A class which encapsulates the functionality of a Cholesky factorization. + /// For a symmetric, positive definite matrix A, the Cholesky factorization + /// is an lower triangular matrix L so that A = L*L'. + /// + /// + /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric + /// or positive definite, the constructor will throw an exception. + /// + public abstract class Cholesky : Cholesky + { + /// + /// Gets the determinant of the matrix for which the Cholesky matrix was computed. + /// + public override float Determinant + { + get + { + var det = 1.0f; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; + } + + return det; + } + } + + /// + /// Gets the log determinant of the matrix for which the Cholesky matrix was computed. + /// + public override float DeterminantLn + { + get + { + var det = 0.0f; + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + det += 2.0f * Convert.ToSingle(Math.Log(CholeskyFactor[j, j])); + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs index 2463a627..e885d196 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -44,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 DenseCholesky : Cholesky + public class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -109,13 +108,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -158,13 +157,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -174,39 +173,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dfactor = (DenseMatrix)CholeskyFactor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override float AddT(float val1, float val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override float LogT(float val1) - { - return (float)Math.Log(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs index 712c4ff2..8ff5bcc7 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -49,16 +48,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class DenseEvd : Evd + public class DenseEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public DenseEvd(DenseMatrix matrix) { @@ -1223,16 +1222,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs index fbe127ac..4b4fb787 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; using Threading; @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class DenseGramSchmidt : GramSchmidt + public class DenseGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an orthogonal matrix @@ -153,13 +152,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -198,41 +197,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs index f4dcebe2..8d199837 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class DenseLU : LU + public class DenseLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -110,13 +109,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment."); } // Copy the contents of input to result. @@ -159,13 +158,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment."); + throw new NotSupportedException("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."); + throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment."); } // Copy the contents of input to result. @@ -186,19 +185,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization Control.LinearAlgebraProvider.LUInverseFactored(result.Data, result.RowCount, Pivots); return result; } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs index 6936dd21..a0fabfc2 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -44,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class DenseQR : QR + public class DenseQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -109,13 +108,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); @@ -154,41 +153,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs index d791c857..6c9208e6 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class DenseSvd : Svd + public class DenseSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -56,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public DenseSvd(DenseMatrix matrix, bool computeVectors) { @@ -117,13 +116,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseMatrix; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } var dresult = result as DenseMatrix; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, input.ColumnCount, dresult.Data); @@ -167,40 +166,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); + throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs new file mode 100644 index 00000000..837d1e5f --- /dev/null +++ b/src/Numerics/LinearAlgebra/Single/Factorization/Evd.cs @@ -0,0 +1,115 @@ +// +// 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.Single.Factorization +{ + using System; + using System.Numerics; + using Generic.Factorization; + + /// + /// Eigenvalues and eigenvectors of a real matrix. + /// + /// + /// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is + /// diagonal and the eigenvector matrix V is orthogonal. + /// I.e. A = V*D*V' and V*VT=I. + /// If A is not symmetric, then the eigenvalue matrix D is block diagonal + /// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, + /// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The + /// columns of V represent the eigenvectors in the sense that A*V = V*D, + /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly + /// conditioned, or even singular, so the validity of the equation + /// A = V*D*Inverse(V) depends upon V.Condition(). + /// + public abstract class Evd : Evd + { + /// + /// Gets the absolute value of determinant of the square matrix for which the EVD was computed. + /// + public override float Determinant + { + get + { + var det = Complex.One; + for (var i = 0; i < VectorEv.Count; i++) + { + det *= VectorEv[i]; + + if (((Numerics.Complex32)VectorEv[i]).AlmostEqual(Numerics.Complex32.Zero)) + { + return 0; + } + } + + return Convert.ToSingle(det.Magnitude); + } + } + + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + var rank = 0; + for (var i = 0; i < VectorEv.Count; i++) + { + if (((Numerics.Complex32)VectorEv[i]).AlmostEqual(Numerics.Complex32.Zero)) + { + continue; + } + + rank++; + } + + return rank; + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < VectorEv.Count; i++) + { + if (VectorEv[i].AlmostEqual(Complex.Zero)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/GramSchmidt.cs b/src/Numerics/LinearAlgebra/Single/Factorization/GramSchmidt.cs new file mode 100644 index 00000000..2a81948e --- /dev/null +++ b/src/Numerics/LinearAlgebra/Single/Factorization/GramSchmidt.cs @@ -0,0 +1,88 @@ +// +// 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.Single.Factorization +{ + using System; + using Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization. + /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix. + /// + /// + /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. + /// + public abstract class GramSchmidt : GramSchmidt + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override float Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) + { + return 0; + } + } + + return Convert.ToSingle(Math.Abs(det)); + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/LU.cs new file mode 100644 index 00000000..36de622d --- /dev/null +++ b/src/Numerics/LinearAlgebra/Single/Factorization/LU.cs @@ -0,0 +1,67 @@ +// +// 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.Single.Factorization +{ + using Generic.Factorization; + + /// + /// 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 : LU + { + /// + /// Gets the determinant of the matrix for which the LU factorization was computed. + /// + public override float Determinant + { + get + { + var det = 1.0f; + for (var j = 0; j < Factors.RowCount; j++) + { + if (Pivots[j] != j) + { + det *= -Factors.At(j, j); + } + else + { + det *= Factors.At(j, j); + } + } + + return det; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/QR.cs new file mode 100644 index 00000000..1b3b6af0 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Single/Factorization/QR.cs @@ -0,0 +1,90 @@ +// +// 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.Single.Factorization +{ + using System; + using Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition. + /// Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m) + /// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix + /// (also called right triangular matrix). + /// + /// + /// The computation of the QR decomposition is done at construction time by Householder transformation. + /// + public abstract class QR : QR + { + /// + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. + /// + public override float Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) + { + return 0; + } + } + + return Convert.ToSingle(Math.Abs(det)); + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public override bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) + { + return false; + } + } + + return true; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs new file mode 100644 index 00000000..61016cfb --- /dev/null +++ b/src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs @@ -0,0 +1,118 @@ +// +// 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.Single.Factorization +{ + using System; + using System.Linq; + using Generic; + using Generic.Factorization; + using Properties; + + /// + /// A class which encapsulates the functionality of the singular value decomposition (SVD). + /// Suppose M is an m-by-n matrix whose entries are real numbers. + /// Then there exists a factorization of the form M = UΣVT where: + /// - U is an m-by-m unitary matrix; + /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; + /// - VT denotes transpose of V, an n-by-n unitary matrix; + /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal + /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined + /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M. + /// + /// + /// The computation of the singular value decomposition is done at construction time. + /// + public abstract class Svd : Svd + { + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public override int Rank + { + get + { + return VectorS.Count(t => !Math.Abs(t).AlmostEqual(0.0f)); + } + } + + /// + /// Gets the two norm of the . + /// + /// The 2-norm of the . + public override float Norm2 + { + get + { + return Math.Abs(VectorS[0]); + } + } + + /// + /// Gets the condition number max(S) / min(S) + /// + /// The condition number. + public override float ConditionNumber + { + get + { + var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; + return Math.Abs(VectorS[0]) / Math.Abs(VectorS[tmp]); + } + } + + /// + /// Gets the determinant of the square matrix for which the SVD was computed. + /// + public override float Determinant + { + get + { + if (MatrixU.RowCount != MatrixVT.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + foreach (var value in VectorS) + { + det *= value; + if (Math.Abs(value).AlmostEqual(0.0f)) + { + return 0; + } + } + + return Convert.ToSingle(Math.Abs(det)); + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs index f2aa0719..6035fa83 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -44,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 class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -220,40 +219,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization result[i] = sum / CholeskyFactor.At(i, i); } } - - #region Simple T Mathematics - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override float AddT(float val1, float val2) - { - return val1 + val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the natural (base e) logarithm of a specified number. - /// - /// A number whose logarithm is to be found - /// Natural (base e) logarithm - protected sealed override float LogT(float val1) - { - return (float)Math.Log(val1); - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs index 8f9aef8a..d1f3c4c2 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization using System; using System.Numerics; using Generic; - using Generic.Factorization; using Numerics; using Properties; @@ -49,16 +48,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// columns of V represent the eigenvectors in the sense that A*V = V*D, /// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly /// conditioned, or even singular, so the validity of the equation - /// A = V*D*Inverse(V) depends upon V.cond(). + /// A = V*D*Inverse(V) depends upon V.Condition(). /// - public class UserEvd : Evd + public class UserEvd : Evd { /// /// Initializes a new instance of the class. This object will compute the /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If EVD algorithm failed to converge with matrix . public UserEvd(Matrix matrix) { @@ -1219,16 +1218,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new ArgumentException(Resources.ArgumentMatrixSymmetric); } } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } } } \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserGramSchmidt.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserGramSchmidt.cs index 44a20708..89453454 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserGramSchmidt.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -42,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization. /// - public class UserGramSchmidt : GramSchmidt + public class UserGramSchmidt : GramSchmidt { /// /// Initializes a new instance of the class. This object creates an orthogonal matrix @@ -69,7 +68,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (var k = 0; k < MatrixQ.ColumnCount; k++) { - var norm = (float)MatrixQ.Column(k).Norm(2); + var norm = MatrixQ.Column(k).Norm(2); if (norm == 0.0) { throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient); @@ -244,30 +243,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs index 4be798b9..2fbf8294 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -43,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public class UserLU : LU + public class UserLU : LU { /// /// Initializes a new instance of the class. This object will compute the @@ -298,19 +297,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization return Solve(inverse); } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs index 5c573594..72a1a947 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization using System; using System.Linq; using Generic; - using Generic.Factorization; using Properties; /// @@ -45,7 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the QR decomposition is done at construction time by Householder transformation. /// - public class UserQR : QR + public class UserQR : QR { /// /// Initializes a new instance of the class. This object will compute the @@ -91,7 +90,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Generate column from initial matrix to work array /// /// Initial matrix - /// The firts row + /// The first row /// The last row /// Column index /// Generated vector @@ -329,29 +328,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization result[i] = inputCopy[i]; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.cs index 8cba9b96..aead860f 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserSvd.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; using Generic; - using Generic.Factorization; using Properties; /// @@ -48,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The computation of the singular value decomposition is done at construction time. /// - public class UserSvd : Svd + public class UserSvd : Svd { /// /// Initializes a new instance of the class. This object will compute the @@ -56,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// /// The matrix to factor. /// Compute the singular U and VT vectors or not. - /// If is null. + /// If is null. /// If SVD algorithm failed to converge with matrix . public UserSvd(Matrix matrix, bool computeVectors) { @@ -478,7 +477,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var shift = 0.0f; if (b != 0.0 || c != 0.0) { - shift = (float) Math.Sqrt((b * b) + c); + shift = (float)Math.Sqrt((b * b) + c); if (b < 0.0) { shift = -shift; @@ -702,14 +701,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization db = z; } - /// dded + /// /// Calculate Norm 2 of the column in matrix starting from row /// /// Source matrix /// The number of rows in /// Column index /// Start row index - /// Norm2 (Euclidean norm) of trhe column + /// Norm2 (Euclidean norm) of the column private static float Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) { float s = 0; @@ -921,30 +920,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization result[j] = value; } } - - #region Simple arithmetic of type T - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Returns the absolute value of a specified number. - /// - /// A number whose absolute is to be found - /// Absolute value - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - - #endregion } } diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b3552989..5f56e271 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -99,8 +99,14 @@ + + + + + + @@ -111,14 +117,26 @@ Code + + + + + + + + + + + + @@ -200,12 +218,18 @@ + + + + + + @@ -242,10 +266,6 @@ - - - - diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index 3240c344..1e8cbf48 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -281,6 +281,9 @@ LinearAlgebra\Complex32\DiagonalMatrix.cs + + LinearAlgebra\Complex32\Factorization\Cholesky.cs + LinearAlgebra\Complex32\Factorization\DenseCholesky.cs @@ -299,6 +302,21 @@ LinearAlgebra\Complex32\Factorization\DenseSvd.cs + + LinearAlgebra\Complex32\Factorization\Evd.cs + + + LinearAlgebra\Complex32\Factorization\GramSchmidt.cs + + + LinearAlgebra\Complex32\Factorization\LU.cs + + + LinearAlgebra\Complex32\Factorization\QR.cs + + + LinearAlgebra\Complex32\Factorization\Svd.cs + LinearAlgebra\Complex32\Factorization\UserCholesky.cs @@ -383,6 +401,9 @@ LinearAlgebra\Complex\DiagonalMatrix.cs + + LinearAlgebra\Complex\Factorization\Cholesky.cs + LinearAlgebra\Complex\Factorization\DenseCholesky.cs @@ -401,6 +422,21 @@ LinearAlgebra\Complex\Factorization\DenseSvd.cs + + LinearAlgebra\Complex\Factorization\Evd.cs + + + LinearAlgebra\Complex\Factorization\GramSchmidt.cs + + + LinearAlgebra\Complex\Factorization\LU.cs + + + LinearAlgebra\Complex\Factorization\QR.cs + + + LinearAlgebra\Complex\Factorization\Svd.cs + LinearAlgebra\Complex\Factorization\UserCholesky.cs @@ -485,6 +521,9 @@ LinearAlgebra\Double\DiagonalMatrix.cs + + LinearAlgebra\Double\Factorization\Cholesky.cs + LinearAlgebra\Double\Factorization\DenseCholesky.cs @@ -503,17 +542,20 @@ LinearAlgebra\Double\Factorization\DenseSvd.cs - - LinearAlgebra\Double\Factorization\SparseCholesky.cs + + LinearAlgebra\Double\Factorization\Evd.cs - - LinearAlgebra\Double\Factorization\SparseLU.cs + + LinearAlgebra\Double\Factorization\GramSchmidt.cs - - LinearAlgebra\Double\Factorization\SparseQR.cs + + LinearAlgebra\Double\Factorization\LU.cs - - LinearAlgebra\Double\Factorization\SparseSvd.cs + + LinearAlgebra\Double\Factorization\QR.cs + + + LinearAlgebra\Double\Factorization\Svd.cs LinearAlgebra\Double\Factorization\UserCholesky.cs @@ -608,6 +650,9 @@ LinearAlgebra\Single\DiagonalMatrix.cs + + LinearAlgebra\Single\Factorization\Cholesky.cs + LinearAlgebra\Single\Factorization\DenseCholesky.cs @@ -626,6 +671,21 @@ LinearAlgebra\Single\Factorization\DenseSvd.cs + + LinearAlgebra\Single\Factorization\Evd.cs + + + LinearAlgebra\Single\Factorization\GramSchmidt.cs + + + LinearAlgebra\Single\Factorization\LU.cs + + + LinearAlgebra\Single\Factorization\QR.cs + + + LinearAlgebra\Single\Factorization\Svd.cs + LinearAlgebra\Single\Factorization\UserCholesky.cs diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs index cba81eae..852e7c2e 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs @@ -55,15 +55,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var I = DenseMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -80,15 +80,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -113,14 +113,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().ConjugateTranspose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().ConjugateTranspose(); for (var i = 0; i < matrix.RowCount; i++) { diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs index 7d2b7bc2..f56d295c 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs @@ -54,15 +54,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var I = UserDefinedMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -79,15 +79,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -112,14 +112,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianUserDefinedMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().ConjugateTranspose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().ConjugateTranspose(); for (var i = 0; i < matrix.RowCount; i++) { diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs index 9adede82..88a12eaf 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs @@ -55,15 +55,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var I = DenseMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -80,15 +80,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -114,14 +114,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().ConjugateTranspose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().ConjugateTranspose(); for (var i = 0; i < matrix.RowCount; i++) { @@ -178,7 +178,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization { var I = DenseMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(1.0, factorEvd.Determinant); + Assert.AreEqual(Numerics.Complex32.One, factorEvd.Determinant); } [Test] diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs index 35f561dd..ac7b0162 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs @@ -54,15 +54,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var I = UserDefinedMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -79,15 +79,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -113,14 +113,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianUserDefinedMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().ConjugateTranspose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().ConjugateTranspose(); for (var i = 0; i < matrix.RowCount; i++) { @@ -177,7 +177,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization { var I = UserDefinedMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(1.0, factorEvd.Determinant); + Assert.AreEqual(Numerics.Complex32.One, factorEvd.Determinant); } [Test] diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs index 156b6d67..29d7d1e6 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs @@ -45,7 +45,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 var value = new Complex32(real, imaginary); var matrix = TestMatrices["Singular3x3"]; var clone = matrix.Clone(); - clone.Multiply(value); + clone = clone.Multiply(value); for (var i = 0; i < matrix.RowCount; i++) { @@ -241,7 +241,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 var matrixB = TestMatrices[mtxB]; var matrix = matrixA.Clone(); - matrix.Add(matrixB); + matrix = matrix.Add(matrixB); for (var i = 0; i < matrix.RowCount; i++) { for (var j = 0; j < matrix.ColumnCount; j++) @@ -341,7 +341,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 var matrixB = TestMatrices[mtxB]; var matrix = matrixA.Clone(); - matrix.Subtract(matrixB); + matrix = matrix.Subtract(matrixB); for (var i = 0; i < matrix.RowCount; i++) { for (var j = 0; j < matrix.ColumnCount; j++) @@ -590,7 +590,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 var matrix = TestMatrices[name]; var copy = matrix.Clone(); - copy.Negate(); + copy = copy.Negate(); for (var i = 0; i < matrix.RowCount; i++) { @@ -794,26 +794,24 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 [Test] public virtual void PointwiseDivideResult() { - foreach (var data in TestMatrices.Values) + var data = TestMatrices["Singular3x3"]; + var other = data.Clone(); + var result = data.Clone(); + data.PointwiseDivide(other, result); + for (var i = 0; i < data.RowCount; i++) { - var other = data.Clone(); - var result = data.Clone(); - data.PointwiseDivide(other, result); - for (var i = 0; i < data.RowCount; i++) + for (var j = 0; j < data.ColumnCount; j++) { - for (var j = 0; j < data.ColumnCount; j++) - { - AssertHelpers.AreEqual(data[i, j] / other[i, j], result[i, j]); - } + AssertHelpers.AreEqual(data[i, j] / other[i, j], result[i, j]); } + } - result = data.PointwiseDivide(other); - for (var i = 0; i < data.RowCount; i++) + result = data.PointwiseDivide(other); + for (var i = 0; i < data.RowCount; i++) + { + for (var j = 0; j < data.ColumnCount; j++) { - for (var j = 0; j < data.ColumnCount; j++) - { - AssertHelpers.AreEqual(data[i, j] / other[i, j], result[i, j]); - } + AssertHelpers.AreEqual(data[i, j] / other[i, j], result[i, j]); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs index 090433e2..09028c39 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs @@ -55,15 +55,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var I = DenseMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -80,15 +80,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -112,14 +112,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().Transpose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().Transpose(); for (var i = 0; i < matrix.RowCount; i++) { diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs index 82b6686a..85ea91fa 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserEvdTests.cs @@ -54,15 +54,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var I = UserDefinedMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -79,15 +79,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -111,14 +111,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().Transpose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().Transpose(); for (var i = 0; i < matrix.RowCount; i++) { diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs index 4e66ed4d..3c79444b 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs @@ -43,7 +43,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { var matrix = TestMatrices["Singular3x3"]; var clone = matrix.Clone(); - clone.Multiply(scalar); + clone = clone.Multiply(scalar); for (var i = 0; i < matrix.RowCount; i++) { @@ -236,7 +236,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double var B = TestMatrices[mtxB]; var matrix = A.Clone(); - matrix.Add(B); + matrix = matrix.Add(B); for (var i = 0; i < matrix.RowCount; i++) { for (var j = 0; j < matrix.ColumnCount; j++) @@ -336,7 +336,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double var B = TestMatrices[mtxB]; var matrix = A.Clone(); - matrix.Subtract(B); + matrix = matrix.Subtract(B); for (var i = 0; i < matrix.RowCount; i++) { for (var j = 0; j < matrix.ColumnCount; j++) @@ -585,7 +585,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double var matrix = TestMatrices[name]; var copy = matrix.Clone(); - copy.Negate(); + copy = copy.Negate(); for (var i = 0; i < matrix.RowCount; i++) { @@ -788,26 +788,24 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public virtual void PointwiseDivideResult() { - foreach (var data in TestMatrices.Values) + var data = TestMatrices["Singular3x3"]; + var other = data.Clone(); + var result = data.Clone(); + data.PointwiseDivide(other, result); + for (var i = 0; i < data.RowCount; i++) { - var other = data.Clone(); - var result = data.Clone(); - data.PointwiseDivide(other, result); - for (var i = 0; i < data.RowCount; i++) + for (var j = 0; j < data.ColumnCount; j++) { - for (var j = 0; j < data.ColumnCount; j++) - { - Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); - } + Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); } + } - result = data.PointwiseDivide(other); - for (var i = 0; i < data.RowCount; i++) + result = data.PointwiseDivide(other); + for (var i = 0; i < data.RowCount; i++) + { + for (var j = 0; j < data.ColumnCount; j++) { - for (var j = 0; j < data.ColumnCount; j++) - { - Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); - } + Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs index 9d3fde4e..76ada795 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs @@ -55,15 +55,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var I = DenseMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -80,15 +80,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -113,14 +113,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().Transpose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().Transpose(); for (var i = 0; i < matrix.RowCount; i++) { @@ -164,7 +164,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization } var factorEvd = matrixA.Evd(); - Assert.AreEqual(factorEvd.Determinant, 0); + AssertHelpers.AlmostEqual(factorEvd.Determinant, 0, 6); Assert.AreEqual(factorEvd.Rank, order - 1); } diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs index e37f6bb8..a007aa9b 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/UserEvdTests.cs @@ -54,15 +54,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var I = UserDefinedMatrix.Identity(order); var factorEvd = I.Evd(); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); - Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(I.RowCount, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); - for (var i = 0; i < factorEvd.EValues().Count; i++) + for (var i = 0; i < factorEvd.EigenValues().Count; i++) { - Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); + Assert.AreEqual(Complex.One, factorEvd.EigenValues()[i]); } } @@ -79,15 +79,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A*V = λ*V - var matrixAv = matrixA * factorEvd.EVectors(); - var matrixLv = factorEvd.EVectors() * factorEvd.D(); + var matrixAv = matrixA * factorEvd.EigenVectors(); + var matrixLv = factorEvd.EigenVectors() * factorEvd.D(); for (var i = 0; i < matrixAv.RowCount; i++) { @@ -112,14 +112,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(order); var factorEvd = matrixA.Evd(); - Assert.AreEqual(order, factorEvd.EVectors().RowCount); - Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); + Assert.AreEqual(order, factorEvd.EigenVectors().RowCount); + Assert.AreEqual(order, factorEvd.EigenVectors().ColumnCount); Assert.AreEqual(order, factorEvd.D().RowCount); Assert.AreEqual(order, factorEvd.D().ColumnCount); // Make sure the A = V*λ*VT - var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().Transpose(); + var matrix = factorEvd.EigenVectors() * factorEvd.D() * factorEvd.EigenVectors().Transpose(); for (var i = 0; i < matrix.RowCount; i++) { diff --git a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs index cee17a32..23bde347 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs @@ -43,7 +43,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single { var matrix = TestMatrices["Singular3x3"]; var clone = matrix.Clone(); - clone.Multiply(scalar); + clone = clone.Multiply(scalar); for (var i = 0; i < matrix.RowCount; i++) { @@ -236,7 +236,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single var B = TestMatrices[mtxB]; var matrix = A.Clone(); - matrix.Add(B); + matrix = matrix.Add(B); for (var i = 0; i < matrix.RowCount; i++) { for (var j = 0; j < matrix.ColumnCount; j++) @@ -336,7 +336,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single var B = TestMatrices[mtxB]; var matrix = A.Clone(); - matrix.Subtract(B); + matrix = matrix.Subtract(B); for (var i = 0; i < matrix.RowCount; i++) { for (var j = 0; j < matrix.ColumnCount; j++) @@ -585,7 +585,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single var matrix = TestMatrices[name]; var copy = matrix.Clone(); - copy.Negate(); + copy = copy.Negate(); for (var i = 0; i < matrix.RowCount; i++) { @@ -788,26 +788,24 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single [Test] public virtual void PointwiseDivideResult() { - foreach (var data in TestMatrices.Values) + var data = TestMatrices["Singular3x3"]; + var other = data.Clone(); + var result = data.Clone(); + data.PointwiseDivide(other, result); + for (var i = 0; i < data.RowCount; i++) { - var other = data.Clone(); - var result = data.Clone(); - data.PointwiseDivide(other, result); - for (var i = 0; i < data.RowCount; i++) + for (var j = 0; j < data.ColumnCount; j++) { - for (var j = 0; j < data.ColumnCount; j++) - { - Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); - } + Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); } + } - result = data.PointwiseDivide(other); - for (var i = 0; i < data.RowCount; i++) + result = data.PointwiseDivide(other); + for (var i = 0; i < data.RowCount; i++) + { + for (var j = 0; j < data.ColumnCount; j++) { - for (var j = 0; j < data.ColumnCount; j++) - { - Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); - } + Assert.AreEqual(data[i, j] / other[i, j], result[i, j]); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.cs index 5fe55245..c44c610c 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.cs @@ -1431,52 +1431,52 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single public virtual void FrobeniusNorm() { var matrix = TestMatrices["Square3x3"]; - AssertHelpers.AlmostEqual(10.7777548f, (float)matrix.FrobeniusNorm(), 7); + AssertHelpers.AlmostEqual(10.7777548f, matrix.FrobeniusNorm(), 7); matrix = TestMatrices["Wide2x3"]; - AssertHelpers.AlmostEqual(4.7947888f, (float)matrix.FrobeniusNorm(), 7); + AssertHelpers.AlmostEqual(4.7947888f, matrix.FrobeniusNorm(), 7); matrix = TestMatrices["Tall3x2"]; - AssertHelpers.AlmostEqual(7.5412200f, (float)matrix.FrobeniusNorm(), 7); + AssertHelpers.AlmostEqual(7.5412200f, matrix.FrobeniusNorm(), 7); } [Test] public virtual void InfinityNorm() { var matrix = TestMatrices["Square3x3"]; - Assert.AreEqual(16.5f, (float)matrix.InfinityNorm()); + AssertHelpers.AlmostEqual(16.5f, matrix.InfinityNorm(), 6); matrix = TestMatrices["Wide2x3"]; - Assert.AreEqual(6.6f, (float)matrix.InfinityNorm()); + AssertHelpers.AlmostEqual(6.6f, matrix.InfinityNorm(), 6); matrix = TestMatrices["Tall3x2"]; - Assert.AreEqual(9.9f, (float)matrix.InfinityNorm()); + AssertHelpers.AlmostEqual(9.9f, matrix.InfinityNorm(), 6); } [Test] public virtual void L1Norm() { var matrix = TestMatrices["Square3x3"]; - Assert.AreEqual(12.1f, (float)matrix.L1Norm()); + Assert.AreEqual(12.1f, matrix.L1Norm()); matrix = TestMatrices["Wide2x3"]; - Assert.AreEqual(5.5f, (float)matrix.L1Norm()); + Assert.AreEqual(5.5f, matrix.L1Norm()); matrix = TestMatrices["Tall3x2"]; - Assert.AreEqual(8.8f, (float)matrix.L1Norm()); + Assert.AreEqual(8.8f, matrix.L1Norm()); } [Test] public virtual void L2Norm() { var matrix = TestMatrices["Square3x3"]; - AssertHelpers.AlmostEqual(10.3913473f, (float)matrix.L2Norm(), 7); + AssertHelpers.AlmostEqual(10.3913473f, matrix.L2Norm(), 7); matrix = TestMatrices["Wide2x3"]; - AssertHelpers.AlmostEqual(4.7540849f, (float)matrix.L2Norm(), 7); + AssertHelpers.AlmostEqual(4.7540849f, matrix.L2Norm(), 7); matrix = TestMatrices["Tall3x2"]; - AssertHelpers.AlmostEqual(7.1827270f, (float)matrix.L2Norm(), 7); + AssertHelpers.AlmostEqual(7.1827270f, matrix.L2Norm(), 7); } } } \ No newline at end of file