diff --git a/src/Numerics/LinearAlgebra/Double/SymmetricDenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SymmetricDenseMatrix.cs new file mode 100644 index 00000000..e34219f5 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/SymmetricDenseMatrix.cs @@ -0,0 +1,820 @@ +// +// 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 +{ + using System; + using Distributions; + using Generic; + + using MathNet.Numerics.LinearAlgebra.Storage; + using MathNet.Numerics.LinearAlgebra.Storage.Indexers.Static; + + using Properties; + using Threading; + + /// + /// A Symmetric Matrix class with dense storage. + /// + /// The underlying storage is a one dimensional array in column-major order. + /// The Upper Triangle is stored(it is equal to the Lower Triangle) + [Serializable] + public class SymmetricDenseMatrix : SymmetricMatrix + { + readonly DenseColumnMajorSymmetricMatrixStorage _storage; + + /// + /// Number of rows. + /// + /// Using this instead of the RowCount property to speed up calculating + /// a matrix index in the data array. + readonly int _rowCount; + + /// + /// Number of columns. + /// + /// Using this instead of the ColumnCount property to speed up calculating + /// a matrix index in the data array. + readonly int _columnCount; + + /// + /// Gets the matrix's data. + /// + /// The matrix's data. + readonly double[] _data; + + internal SymmetricDenseMatrix(DenseColumnMajorSymmetricMatrixStorage storage) + : base(storage) + { + _storage = storage; + _rowCount = _storage.RowCount; + _columnCount = _storage.ColumnCount; + _data = _storage.Data; + } + + /// + /// Initializes a new instance of the class. This matrix is square with a given size. + /// + /// The order of the matrix. + /// + /// If is less than one. + /// + public SymmetricDenseMatrix(int order) + : this(new DenseColumnMajorSymmetricMatrixStorage(order)) + { + } + + /// + /// Initializes a new instance of the class with all entries set to a particular value. + /// + /// + /// The order of the matrix. + /// + /// The value which we assign to each element of the matrix. + /// Forcing user to input (int, int, double) because only asking (int, double) would + /// create a signature easily confused with (int, int) which is already used + public SymmetricDenseMatrix(int order, double value) + : this(order) + { + for (var i = 0; i < Data.Length; i++) + { + Data[i] = value; + } + } + + /// + /// Initializes a new instance of the class from a one dimensional array. This constructor + /// will reference the one dimensional array and not copy it. + /// + /// The size of the square matrix. + /// + /// The one dimensional array to create this matrix from. Column-major and row-major order is identical on a symmetric matrix: http://en.wikipedia.org/wiki/Row-major_order + /// + /// + /// If does not represent a packed array. + /// + public SymmetricDenseMatrix(int order, double[] array) + : this(new DenseColumnMajorSymmetricMatrixStorage(order, array)) + { + } + + /// + /// Initializes a new instance of the class from a 2D array. This constructor + /// will allocate a completely new memory block for storing the symmetric dense matrix. + /// + /// The 2D array to create this matrix from. + /// + /// If is not a square array. + /// + /// + /// If is not a symmetric array. + /// + public SymmetricDenseMatrix(double[,] array) + : this(array.GetLength(0)) + { + if (!CheckIfSymmetric(array)) + { + throw new ArgumentException(Resources.ArgumentMatrixSymmetric); + } + + var indexer = new PackedStorageIndexerUpper(Order); + for (var row = 0; row < Order; row++) + { + for (var column = row; column < Order; column++) + { + Data[indexer.Of(row, column)] = array[row, column]; + } + } + } + + /// + /// Initializes a new instance of the class, copying + /// the values from the given matrix. Matrix must be Symmetric. + /// + /// The matrix to copy. + /// + /// If is not a square matrix. + /// + /// + /// If is not a symmetric matrix. + /// + public SymmetricDenseMatrix(Matrix matrix) + : this(matrix.RowCount) + { + var symmetricMatrix = matrix as SymmetricDenseMatrix; + + if (!matrix.IsSymmetric) + { + throw new ArgumentException(Resources.ArgumentMatrixSymmetric); + } + + if (symmetricMatrix == null) + { + var indexer = new PackedStorageIndexerUpper(Order); + for (var row = 0; row < Order; row++) + { + for (var column = row; column < Order; column++) + { + Data[indexer.Of(row, column)] = matrix[row, column]; + } + } + } + else + { + matrix.CopyTo(this); + } + } + + /// + /// Gets the matrix's data in array format. + /// + /// The matrix's raw data. + public double[] Data + { + get; + private set; + } + + /// + /// Creates a SymmetricDenseMatrix for the given number of rows and columns. + /// If rows and columns are not equal, returns a DenseMatrix instead. + /// + /// + /// The number of rows. + /// + /// + /// The number of columns. + /// + /// True if all fields must be mutable (e.g. not a diagonal matrix). + /// + /// A DenseMatrix or SymmetricDenseMatrix with the given dimensions. + /// + /// /// + /// If is not equal to . + /// Symmetric arrays are always square + /// + public override Matrix CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false) + { + if (numberOfRows != numberOfColumns || fullyMutable) + { + return new DenseMatrix(numberOfRows, numberOfColumns); + } + + return new SymmetricDenseMatrix(numberOfRows, numberOfColumns); + } + + /// + /// Creates a with a the given dimension. + /// + /// The size of the vector. + /// True if all fields must be mutable. + /// + /// A with the given dimension. + /// + public override Vector CreateVector(int size, bool fullyMutable = false) + { + return new DenseVector(size); + } + + #region Static constructors for special matrices. + + /// + /// Initializes a square with all zero's except for ones on the diagonal. + /// + /// the size of the square matrix. + /// A symmetric dense identity matrix. + /// + /// If is less than one. + /// + public static SymmetricDenseMatrix Identity(int order) + { + var m = new SymmetricDenseMatrix(order); + for (var i = 0; i < order; i++) + { + m.At(i, i, 1.0); + } + + return m; + } + + #endregion + + /// + /// Adds another matrix to this matrix. + /// + /// The matrix to add to this matrix. + /// The matrix to store the result of add + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + protected override void DoAdd(Matrix other, Matrix result) + { + var denseOther = other as SymmetricDenseMatrix; + var denseResult = result as SymmetricDenseMatrix; + if (denseOther == null || denseResult == null) + { + base.DoAdd(other, result); + } + else + { + Control.LinearAlgebraProvider.AddArrays(Data, denseOther.Data, denseResult.Data); + } + } + + /// + /// Subtracts another matrix from this matrix. + /// + /// The matrix to subtract. + /// The matrix to store the result of the subtraction. + protected override void DoSubtract(Matrix other, Matrix result) + { + var denseOther = other as SymmetricDenseMatrix; + var denseResult = result as SymmetricDenseMatrix; + if (denseOther == null || denseResult == null) + { + base.DoSubtract(other, result); + } + else + { + Control.LinearAlgebraProvider.SubtractArrays(Data, denseOther.Data, denseResult.Data); + } + } + + /// + /// Multiplies each element of the matrix by a scalar and places results into the result matrix. + /// + /// The scalar to multiply the matrix with. + /// The matrix to store the result of the multiplication. + protected override void DoMultiply(double scalar, Matrix result) + { + var denseResult = result as SymmetricDenseMatrix; + if (denseResult == null) + { + base.DoMultiply(scalar, result); + } + else + { + Control.LinearAlgebraProvider.ScaleArray(scalar, Data, denseResult.Data); + } + } + + /// + /// Multiplies this matrix with a vector and places the results into the result vector. + /// + /// The vector to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Vector rightSide, Vector result) + { + var denseRight = rightSide as DenseVector; + var denseResult = result as DenseVector; + + if (denseRight == null || denseResult == null) + { + base.DoMultiply(rightSide, result); + } + else + { + // TODO: Change this when symmetric methods are implemented in the Linear Algebra Providers. + base.DoMultiply(rightSide, result); + } + } + + /// + /// Multiplies this matrix with another matrix and places the results into the result matrix. + /// + /// The matrix to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Matrix other, Matrix result) + { + var denseOther = other as SymmetricDenseMatrix; + var denseResult = result as SymmetricDenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoMultiply(other, result); + } + else + { + // TODO: Change this when symmetric methods are implemented in the Linear Algebra Providers. + base.DoMultiply(other, result); + } + } + + /// + /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix. + /// + /// The matrix to multiply with. + /// The result of the multiplication. + protected override void DoTransposeAndMultiply(Matrix other, Matrix result) + { + var denseOther = other as SymmetricDenseMatrix; + var denseResult = result as SymmetricDenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoTransposeAndMultiply(other, result); + } + else + { + // TODO: Change this when symmetric methods are implemented in the Linear Algebra Providers. + base.DoTransposeAndMultiply(other, result); + } + } + + /// + /// Negate each element of this matrix and place the results into the result matrix. + /// + /// The result of the negation. + protected override void DoNegate(Matrix result) + { + var denseResult = result as SymmetricDenseMatrix; + + if (denseResult == null) + { + base.DoNegate(result); + } + else + { + Control.LinearAlgebraProvider.ScaleArray(-1, Data, denseResult.Data); + } + } + + /// + /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix. + /// + /// The matrix to pointwise multiply with this one. + /// The matrix to store the result of the pointwise multiplication. + protected override void DoPointwiseMultiply(Matrix other, Matrix result) + { + var denseOther = other as SymmetricDenseMatrix; + var denseResult = result as SymmetricDenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoPointwiseMultiply(other, result); + } + else + { + Control.LinearAlgebraProvider.PointWiseMultiplyArrays(Data, denseOther.Data, denseResult.Data); + } + } + + /// + /// Pointwise divide this matrix by another matrix and stores the result into the result matrix. + /// + /// The matrix to pointwise divide this one by. + /// The matrix to store the result of the pointwise division. + protected override void DoPointwiseDivide(Matrix other, Matrix result) + { + var denseOther = other as SymmetricDenseMatrix; + var denseResult = result as SymmetricDenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoPointwiseDivide(other, result); + } + else + { + Control.LinearAlgebraProvider.PointWiseDivideArrays(Data, denseOther.Data, denseResult.Data); + } + } + + /// + /// Returns a new matrix containing the lower triangle of this matrix. + /// + /// The lower triangle of this matrix. + public override Matrix LowerTriangle() + { + var ret = new DenseMatrix(Order); + for (var row = 0; row < Order; row++) + { + for (var column = 0; column <= row; column++) + { + ret[row, column] = At(row, column); + } + } + + return ret; + } + + /// + /// Returns a new matrix containing the lower triangle of this matrix. The new matrix + /// does not contain the diagonal elements of this matrix. + /// + /// The lower triangle of this matrix. + public override Matrix StrictlyLowerTriangle() + { + var ret = new DenseMatrix(Order); + for (var row = 0; row < Order; row++) + { + for (var column = 0; column < row; column++) + { + ret[row, column] = At(row, column); + } + } + + return ret; + } + + /// + /// Returns a new matrix containing the upper triangle of this matrix. + /// + /// The upper triangle of this matrix. + public override Matrix UpperTriangle() + { + var ret = new DenseMatrix(Order); + for (var row = 0; row < Order; row++) + { + for (var column = row; column < Order; column++) + { + ret[row, column] = At(row, column); + } + } + + return ret; + } + + /// + /// Returns a new matrix containing the upper triangle of this matrix. The new matrix + /// does not contain the diagonal elements of this matrix. + /// + /// The upper triangle of this matrix. + public override Matrix StrictlyUpperTriangle() + { + var ret = new DenseMatrix(Order); + for (var row = 0; row < Order; row++) + { + for (var column = row + 1; column < Order; column++) + { + ret[row, column] = At(row, column); + } + } + + return ret; + } + + /// + /// Computes the modulus for each element of the matrix. + /// + /// The divisor to use. + /// Matrix to store the results in. + protected override void DoModulus(double divisor, Matrix result) + { + var denseResult = result as SymmetricDenseMatrix; + + if (denseResult == null) + { + base.DoModulus(divisor, result); + } + else + { + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } + + CommonParallel.For( + 0, + Data.Length, + index => denseResult.Data[index] %= divisor); + } + } + + /// + /// Computes the trace of this matrix. + /// + /// The trace of this matrix + /// If the matrix is not square + public override double Trace() + { + if (RowCount != ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var sum = 0.0; + for (var i = 0; i < RowCount; i++) + { + sum += At(i, i); + } + + return sum; + } + + /// + /// Populates a symmetric matrix with random elements. + /// + /// The symmetric matrix to populate. + /// Continuous Random Distribution to generate elements from. + protected override void DoRandom(Matrix matrix, IContinuousDistribution distribution) + { + var denseMatrix = matrix as SymmetricDenseMatrix; + + if (denseMatrix == null) + { + base.DoRandom(matrix, distribution); + } + else + { + for (var i = 0; i < denseMatrix.Data.Length; i++) + { + denseMatrix.Data[i] = distribution.Sample(); + } + } + } + + /// + /// Populates a symmetric matrix with random elements. + /// + /// The symmetric matrix to populate. + /// Continuous Random Distribution to generate elements from. + protected override void DoRandom(Matrix matrix, IDiscreteDistribution distribution) + { + var denseMatrix = matrix as SymmetricDenseMatrix; + + if (denseMatrix == null) + { + base.DoRandom(matrix, distribution); + } + else + { + for (var i = 0; i < denseMatrix.Data.Length; i++) + { + denseMatrix.Data[i] = distribution.Sample(); + } + } + } + + /// + /// Adds two matrices together and returns the results. + /// + /// This operator will allocate new memory for the result. It will + /// choose the representation of either or depending on which + /// is denser. + /// The left matrix to add. + /// The right matrix to add. + /// The result of the addition. + /// If and don't have the same dimensions. + /// If or is . + public static SymmetricDenseMatrix operator +(SymmetricDenseMatrix leftSide, SymmetricDenseMatrix rightSide) + { + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + if (leftSide == null) + { + throw new ArgumentNullException("leftSide"); + } + + if (leftSide.RowCount != rightSide.RowCount) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); + } + + return (SymmetricDenseMatrix)leftSide.Add(rightSide); + } + + /// + /// Returns a Matrix containing the same values of . + /// + /// The matrix to get the values from. + /// A matrix containing a the same values as . + /// If is . + public static SymmetricDenseMatrix operator +(SymmetricDenseMatrix rightSide) + { + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + return (SymmetricDenseMatrix)rightSide.Clone(); + } + + /// + /// Subtracts two matrices together and returns the results. + /// + /// This operator will allocate new memory for the result. It will + /// choose the representation of either or depending on which + /// is denser. + /// The left matrix to subtract. + /// The right matrix to subtract. + /// The result of the addition. + /// If and don't have the same dimensions. + /// If or is . + public static SymmetricDenseMatrix operator -(SymmetricDenseMatrix leftSide, SymmetricDenseMatrix rightSide) + { + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + if (leftSide == null) + { + throw new ArgumentNullException("leftSide"); + } + + if (leftSide.RowCount != rightSide.RowCount) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); + } + + return (SymmetricDenseMatrix)leftSide.Subtract(rightSide); + } + + /// + /// Negates each element of the matrix. + /// + /// The matrix to negate. + /// A matrix containing the negated values. + /// If is . + public static SymmetricDenseMatrix operator -(SymmetricDenseMatrix rightSide) + { + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + return (SymmetricDenseMatrix)rightSide.Negate(); + } + + /// + /// Multiplies a Matrix by a constant and returns the result. + /// + /// The matrix to multiply. + /// The constant to multiply the matrix by. + /// The result of the multiplication. + /// If is . + public static SymmetricDenseMatrix operator *(SymmetricDenseMatrix leftSide, double rightSide) + { + if (leftSide == null) + { + throw new ArgumentNullException("leftSide"); + } + + return (SymmetricDenseMatrix)leftSide.Multiply(rightSide); + } + + /// + /// Multiplies a Matrix by a constant and returns the result. + /// + /// The matrix to multiply. + /// The constant to multiply the matrix by. + /// The result of the multiplication. + /// If is . + public static SymmetricDenseMatrix operator *(double leftSide, SymmetricDenseMatrix rightSide) + { + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + return (SymmetricDenseMatrix)rightSide.Multiply(leftSide); + } + + /// + /// Multiplies two matrices. + /// + /// This operator will allocate new memory for the result. It will + /// choose the representation of either or depending on which + /// is denser. + /// The left matrix to multiply. + /// The right matrix to multiply. + /// The result of multiplication. + /// If or is . + /// If the dimensions of or don't conform. + public static SymmetricDenseMatrix operator *(SymmetricDenseMatrix leftSide, SymmetricDenseMatrix rightSide) + { + if (leftSide == null) + { + throw new ArgumentNullException("leftSide"); + } + + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + if (leftSide.ColumnCount != rightSide.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + return (SymmetricDenseMatrix)leftSide.Multiply(rightSide); + } + + /// + /// Multiplies a Matrix and a Vector. + /// + /// The matrix to multiply. + /// The vector to multiply. + /// The result of multiplication. + /// If or is . + public static DenseVector operator *(SymmetricDenseMatrix leftSide, DenseVector rightSide) + { + if (leftSide == null) + { + throw new ArgumentNullException("leftSide"); + } + + return (DenseVector)leftSide.Multiply(rightSide); + } + + /// + /// Multiplies a Vector and a Matrix. + /// + /// The vector to multiply. + /// The matrix to multiply. + /// The result of multiplication. + /// If or is . + public static DenseVector operator *(DenseVector leftSide, SymmetricDenseMatrix rightSide) + { + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + return (DenseVector)rightSide.LeftMultiply(leftSide); + } + + /// + /// Multiplies a Matrix by a constant and returns the result. + /// + /// The matrix to multiply. + /// The constant to multiply the matrix by. + /// The result of the multiplication. + /// If is . + public static SymmetricDenseMatrix operator %(SymmetricDenseMatrix leftSide, double rightSide) + { + if (leftSide == null) + { + throw new ArgumentNullException("leftSide"); + } + + return (SymmetricDenseMatrix)leftSide.Modulus(rightSide); + } + } +} \ No newline at end of file diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorSymmetricMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorSymmetricMatrixStorage.cs index a2b52833..3d0ee7f7 100644 --- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorSymmetricMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorSymmetricMatrixStorage.cs @@ -54,7 +54,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage /// Not range-checked. public override T At(int row, int column) { - return Data[Indexer.Of(row, column)]; + var r = Math.Min(row, column); + var c = Math.Max(row, column); + return Data[Indexer.Of(r, c)]; } /// @@ -66,6 +68,11 @@ namespace MathNet.Numerics.LinearAlgebra.Storage /// WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks. public override void At(int row, int column, T value) { + if (row > column) + { + throw new IndexOutOfRangeException("Setting an element in the strictly lower triangle of a symmetric matrix is disabled to avoid errors"); + } + Data[Indexer.Of(row, column)] = value; } diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 531cbd69..035cb821 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -141,6 +141,7 @@ +