From fb943d5522bfaba6b8489e6e6be3e51de6ceefbd Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Thu, 21 Oct 2010 21:08:36 +0800 Subject: [PATCH] vector: cleaned up the vector class a little more. now starting the matrix class and the factorization classes --- .../LinearAlgebra/Double/DenseMatrix.cs | 558 +++++++++--------- src/Numerics/LinearAlgebra/Generic/Common.cs | 64 +- .../Generic/Factorization/Svd.cs | 5 +- .../Generic/Matrix.Arithmetic.cs | 311 +++++----- src/Numerics/LinearAlgebra/Generic/Matrix.cs | 186 +----- src/Numerics/LinearAlgebra/Generic/Vector.cs | 111 +--- src/Numerics/Numerics.csproj | 18 +- 7 files changed, 535 insertions(+), 718 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index ad8acd72..7d9792bd 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -35,7 +35,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// A Matrix class with dense storage. The underlying storage is a one dimensional array in column-major order. /// - public class DenseMatrix : Matrix + public class DenseMatrix : Matrix { /// /// Initializes a new instance of the class. This matrix is square with a given size. @@ -242,7 +242,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override double FrobeniusNorm() { var transpose = (DenseMatrix)Transpose(); - var aat = this * transpose; + var aat = (DenseMatrix) (this * transpose); var norm = 0.0; for (var i = 0; i < RowCount; i++) @@ -276,272 +276,45 @@ namespace MathNet.Numerics.LinearAlgebra.Double #region Elementary operations /// - /// Adds another matrix to this matrix. The result will be written into this matrix. + /// Adds another matrix to this matrix. /// /// The matrix to add to this matrix. - /// If the other matrix is . + /// The matrix to store the result of add + /// If the other matrix is . /// If the two matrices don't have the same dimensions. - public override void Add(Matrix other) + protected override void DoAdd(Matrix other, Matrix result) { - var m = other as DenseMatrix; - if (m == null) + var denseOther = other as DenseMatrix; + var denseResult = result as DenseMatrix; + if (denseOther == null || denseResult == null) { - base.Add(other); + base.DoAdd(other, result); } else { - Add(m); + Control.LinearAlgebraProvider.AddArrays(Data, denseOther.Data, denseResult.Data); } } /// - /// Adds another to this matrix. The result will be written into this matrix. - /// - /// The to add to this matrix. - /// If the other matrix is . - /// If the two matrices don't have the same dimensions. - public void Add(DenseMatrix other) - { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (other.RowCount != RowCount || other.ColumnCount != ColumnCount) - { - throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); - } - - Control.LinearAlgebraProvider.AddArrays(Data, other.Data, Data); - } - - /// - /// Subtracts another matrix from this matrix. The result will be written into this matrix. + /// Subtracts another matrix from this matrix. /// /// The matrix to subtract. - /// If the other matrix is . - /// If the two matrices don't have the same dimensions. - public override void Subtract(Matrix other) - { - var m = other as DenseMatrix; - if (m == null) - { - base.Subtract(other); - } - else - { - Subtract(m); - } - } - - /// - /// Subtracts another from this matrix. The result will be written into this matrix. - /// - /// The to subtract. - /// If the other matrix is . - /// If the two matrices don't have the same dimensions. - public void Subtract(DenseMatrix other) - { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (other.RowCount != RowCount || other.ColumnCount != ColumnCount) - { - throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); - } - - Control.LinearAlgebraProvider.SubtractArrays(Data, other.Data, Data); - } - - /// - /// Multiplies each element of this matrix with a scalar. - /// - /// The scalar to multiply with. - public override void Multiply(double scalar) - { - Control.LinearAlgebraProvider.ScaleArray(scalar, Data); - } - - /// - /// Multiplies this dense matrix with another dense matrix and places the results into the result dense matrix. - /// - /// The matrix to multiply with. - /// The result of the multiplication. - /// If the other matrix is . - /// If the result matrix is . - /// If this.Columns != other.Rows. - /// If the result matrix's dimensions are not the this.Rows x other.Columns. - public override void Multiply(Matrix other, Matrix result) + /// The matrix to store the result of the subtraction. + protected override void DoSubtract(Matrix other, Matrix result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (result == null) + var denseOther = other as DenseMatrix; + var denseResult = result as DenseMatrix; + if (denseOther == null || denseResult == null) { - throw new ArgumentNullException("result"); - } - - if (ColumnCount != other.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - if (result.RowCount != RowCount || result.ColumnCount != other.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - var m = other as DenseMatrix; - var r = result as DenseMatrix; - - if (m == null || r == null) - { - base.Multiply(other, result); + base.DoSubtract(other, result); } else { - Control.LinearAlgebraProvider.MatrixMultiply( - Data, - RowCount, - ColumnCount, - m.Data, - m.RowCount, - m.ColumnCount, - r.Data); - } - } - - /// - /// Multiplies this matrix with another matrix and returns the result. - /// - /// The matrix to multiply with. - /// If this.Columns != other.Rows. - /// If the other matrix is . - /// The result of multiplication. - public override Matrix Multiply(Matrix other) - { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (ColumnCount != other.RowCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - var m = other as DenseMatrix; - if (m == null) - { - return base.Multiply(other); - } - - var result = (DenseMatrix)CreateMatrix(RowCount, other.ColumnCount); - Multiply(other, result); - return result; - } - - /// - /// Multiplies this dense matrix with transpose of another dense matrix and places the results into the result dense matrix. - /// - /// The matrix to multiply with. - /// The result of the multiplication. - /// If the other matrix is . - /// If the result matrix is . - /// If this.Columns != other.Rows. - /// If the result matrix's dimensions are not the this.Rows x other.Columns. - public override void TransposeAndMultiply(Matrix other, Matrix result) - { - var otherDense = other as DenseMatrix; - var resultDense = result as DenseMatrix; - - if (otherDense == null || resultDense == null) - { - base.TransposeAndMultiply(other, result); - return; - } - - if (ColumnCount != otherDense.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - if ((resultDense.RowCount != RowCount) || (resultDense.ColumnCount != otherDense.RowCount)) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); - } - - Control.LinearAlgebraProvider.MatrixMultiplyWithUpdate( - Algorithms.LinearAlgebra.Transpose.DontTranspose, - Algorithms.LinearAlgebra.Transpose.Transpose, - 1.0, - Data, - RowCount, - ColumnCount, - otherDense.Data, - otherDense.RowCount, - otherDense.ColumnCount, - 1.0, - resultDense.Data); - } - - /// - /// Multiplies this matrix with transpose of another matrix and returns the result. - /// - /// The matrix to multiply with. - /// If this.Columns != other.Rows. - /// If the other matrix is . - /// The result of multiplication. - public override Matrix TransposeAndMultiply(Matrix other) - { - var otherDense = other as DenseMatrix; - if (otherDense == null) - { - return base.TransposeAndMultiply(other); - } - - if (ColumnCount != otherDense.ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixDimensions); + Control.LinearAlgebraProvider.SubtractArrays(Data, denseOther.Data, denseResult.Data); } - - var result = (DenseMatrix)CreateMatrix(RowCount, other.RowCount); - TransposeAndMultiply(other, result); - return result; - } - - /// - /// Multiplies two dense matrices. - /// - /// 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 DenseMatrix operator *(DenseMatrix leftSide, DenseMatrix 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 (DenseMatrix)leftSide.Multiply(rightSide); } - + #endregion #region Static constructors for special matrices. @@ -567,14 +340,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double #endregion - /// - /// Negates each element of this matrix. - /// - public override void Negate() - { - Multiply(-1); - } - /// /// Generates matrix with random elements. /// @@ -651,61 +416,280 @@ namespace MathNet.Numerics.LinearAlgebra.Double return matrix; } - #region Simple arithmetic of type T /// - /// Add two values T+T + /// Returns the conjugate transpose of this matrix. + /// + /// The conjugate transpose of this matrix. + public override Matrix ConjugateTranspose() + { + return Transpose(); + } + + /* Control.LinearAlgebraProvider.MatrixMultiplyWithUpdate( + Algorithms.LinearAlgebra.Transpose.DontTranspose, + Algorithms.LinearAlgebra.Transpose.Transpose, + 1.0, + Data, + RowCount, + ColumnCount, + otherDense.Data, + otherDense.RowCount, + otherDense.ColumnCount, + 1.0, + resultDense.Data); + */ + + /// + /// Multiplies each element of the matrix by a scalar and places results into the result matrix. /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override double AddT(double val1, double val2) + /// The scalar to multiply the matrix with. + /// The matrix to store the result of the multiplication. + protected override void DoMultiply(double scalar, Matrix result) { - return val1 + val2; + var denseResult = result as DenseMatrix; + if (denseResult == null) + { + base.DoMultiply(scalar, result); + } + else + { + Control.LinearAlgebraProvider.ScaleArray(scalar, Data); + } } /// - /// Subtract two values T-T + /// Multiplies this matrix with a vector and places the results into the result vector. /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override double SubtractT(double val1, double val2) + /// The vector to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Vector rightSide, Vector result) { - return val1 - val2; + CommonParallel.For( + 0, + RowCount, + i => + { + var s = 0.0; + for (var j = 0; j != ColumnCount; j++) + { + s += At(i, j) * rightSide[j]; + } + + result[i] = s; + }); } /// - /// Multiply two values T*T + /// Left multiply a matrix with a vector ( = vector * matrix ) and place the result in the result vector. /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) + /// The vector to multiply with. + /// The result of the multiplication. + protected override void DoLeftMultiply(Vector leftSide, Vector result) { - return val1 * val2; + CommonParallel.For( + 0, + RowCount, + j => + { + var s = 0.0; + for (var i = 0; i != leftSide.Count; i++) + { + s += leftSide[i] * At(i, j); + } + + result[j] = s; + }); } /// - /// Divide two values T/T + /// Multiplies this matrix with another matrix and places the results into the result matrix. /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override double DivideT(double val1, double val2) + /// The matrix to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Matrix other, Matrix result) { - return val1 / val2; + var denseOther = other as DenseMatrix; + var denseResult = result as DenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoMultiply(other, result); + } + else + { + CommonParallel.For( + 0, + RowCount, + j => + { + for (var i = 0; i != other.ColumnCount; i++) + { + var s = 0.0; + for (var l = 0; l < ColumnCount; l++) + { + s += Data[(j * RowCount) + l] * denseOther.Data[(i * RowCount) + l]; + } + + result.At(j, i, s); + } + }); + + CommonParallel.For( + 0, + RowCount, + j => + { + for (var i = 0; i < RowCount; i++) + { + var s = 0.0; + for (var l = 0; l < ColumnCount; l++) + { + s += Data[(j * RowCount) + l] * denseOther.Data[(l * RowCount) + j]; + } + + denseResult.Data[(j * RowCount) + i] *= s; + } + }); + } + } + + /// + /// 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 DenseMatrix; + var denseResult = result as DenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoTransposeAndMultiply(other, result); + } + else + { + CommonParallel.For( + 0, + RowCount, + j => + { + for (var i = 0; i < RowCount; i++) + { + var s = 0.0; + for (var l = 0; l < ColumnCount; l++) + { + s += Data[(j * RowCount) + l] * denseOther.Data[(l * RowCount) + j]; + } + + denseResult.Data[(j * RowCount) + i] *= s; + } + }); + } } /// - /// Take absolute value + /// Negate each element of this matrix and place the results into the result matrix. /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(double val1) + /// The result of the negation. + protected override void DoNegate(Matrix result) { - return Math.Abs(val1); + var denseResult = result as DenseMatrix; + + if (denseResult == null) + { + base.DoNegate(result); + } + else + { + CommonParallel.For( + 0, + RowCount, + i => + { + for (var j = 0; j != ColumnCount; j++) + { + var index = (j * RowCount) + i; + denseResult.Data[index] =- Data[index]; + } + }); + } + } + + /// + /// 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 DenseMatrix; + var denseResult = result as DenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoPointwiseMultiply(other, result); + } + else + { + CommonParallel.For( + 0, + ColumnCount, + j => + { + for (var i = 0; i < RowCount; i++) + { + var index = (j * RowCount) + i; + denseResult.Data[index] = Data[index] * denseOther.Data[index]; + + } + }); + } + } + + /// + /// 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 DenseMatrix; + var denseResult = result as DenseMatrix; + + if (denseOther == null || denseResult == null) + { + base.DoPointwiseDivide(other, result); + } + else + { + CommonParallel.For( + 0, + ColumnCount, + j => + { + for (var i = 0; i < RowCount; i++) + { + var index = (j * RowCount) + i; + denseResult.Data[index] = Data[index] / denseOther.Data[index]; + } + }); + } + } + + /// + /// 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); + } + + return CommonParallel.Aggregate(0, RowCount, i => Data[(i * RowCount) + i]); } - #endregion - } } diff --git a/src/Numerics/LinearAlgebra/Generic/Common.cs b/src/Numerics/LinearAlgebra/Generic/Common.cs index ae443530..124b1be6 100644 --- a/src/Numerics/LinearAlgebra/Generic/Common.cs +++ b/src/Numerics/LinearAlgebra/Generic/Common.cs @@ -1,10 +1,36 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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.Generic { + using System; + + /// + /// A setup functions to help simplify the generic code. + /// internal static class Common { /// @@ -17,5 +43,35 @@ namespace MathNet.Numerics.LinearAlgebra.Generic { return Math.Max(a, b); } + + /// + /// Sets the value of 1.0 for type T. + /// + /// The type to return the value of 1.0 of. + /// The value of 1.0 for type T. + public static T SetOne() + { + if (typeof(T) == typeof(System.Numerics.Complex)) + { + return (T)(object)System.Numerics.Complex.One; + } + + if (typeof(T) == typeof(Numerics.Complex32)) + { + return (T)(object)Numerics.Complex32.One; + } + + if (typeof(T) == typeof(double)) + { + return (T)(object)1.0; + } + + if (typeof(T) == typeof(float)) + { + return (T)(object)1.0f; + } + + throw new NotSupportedException(); + } } } diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs index 9d3774de..1ba48da3 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/Svd.cs @@ -162,11 +162,12 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// Gets the two norm of the . /// /// The 2-norm of the . - public virtual double Norm2 + public virtual T Norm2 { get { - return AbsoluteT(VectorS[0]); + throw new NotImplementedException(); + //return AbsoluteT(VectorS[0]); } } diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs index cf75523d..e819a872 100644 --- a/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs @@ -40,12 +40,23 @@ namespace MathNet.Numerics.LinearAlgebra.Generic public abstract partial class Matrix { /// - /// Adds another matrix to this matrix. The result will be written into this matrix. + /// The value of 1.0. + /// + private static readonly T One = Common.SetOne(); + + /// + /// The value of 0.0. + /// + private static readonly T Zero = default(T); + + /// + /// Adds another matrix to this matrix. /// /// The matrix to add to this matrix. - /// If the other matrix is . + /// The result of the addition. + /// If the other matrix is . /// If the two matrices don't have the same dimensions. - public virtual void Add(Matrix other) + public virtual Matrix Add(Matrix other) { if (other == null) { @@ -57,25 +68,55 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); } - CommonParallel.For( - 0, - RowCount, - i => - { - for (var j = 0; j < ColumnCount; j++) - { - At(i, j, AddT(At(i, j), other.At(i, j))); - } - }); + var result = CreateMatrix(RowCount, ColumnCount); + Add(other, result); + return result; } /// - /// Subtracts another matrix from this matrix. The result will be written into this matrix. + /// Adds another matrix to this matrix. + /// + /// The matrix to add to this matrix. + /// The matrix to store the result of the addition. + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + public virtual void Add(Matrix other, Matrix result) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (other.RowCount != RowCount || other.ColumnCount != ColumnCount) + { + throw new ArgumentOutOfRangeException("other", Resources.ArgumentMatrixDimensions); + } + + if (result.RowCount != RowCount || result.ColumnCount != ColumnCount) + { + throw new ArgumentOutOfRangeException("result", Resources.ArgumentMatrixDimensions); + } + + DoAdd(other, result); + } + + /// + /// Adds another matrix to this matrix. + /// + /// The matrix to add to this matrix. + /// The matrix to store the result of the addition. + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + protected abstract void DoAdd(Matrix other, Matrix result); + + /// + /// Subtracts another matrix from this matrix. /// /// The matrix to subtract. - /// If the other matrix is . + /// The result of the subtraction. + /// If the other matrix is . /// If the two matrices don't have the same dimensions. - public virtual void Subtract(Matrix other) + public virtual Matrix Subtract(Matrix other) { if (other == null) { @@ -87,46 +128,67 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); } - CommonParallel.For( - 0, - RowCount, - i => - { - for (var j = 0; j < ColumnCount; j++) - { - At(i, j, SubtractT(At(i, j), other.At(i, j))); - } - }); + var result = CreateMatrix(RowCount, ColumnCount); + DoSubtract(other, result); + return result; + } + + /// + /// Subtracts another matrix from this matrix. + /// + /// The matrix to subtract. + /// The matrix to store the result of the subtraction. + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + public virtual void Subtract(Matrix other, Matrix result) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (other.RowCount != RowCount || other.ColumnCount != ColumnCount) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); + } + + DoSubtract(other, result); } + /// + /// Subtracts another matrix from this matrix. + /// + /// The matrix to subtract. + /// The matrix to store the result of the subtraction. + protected abstract void DoSubtract(Matrix other, Matrix result); + /// /// Multiplies each element of this matrix with a scalar. /// /// The scalar to multiply with. - public virtual void Multiply(T scalar) + /// The result of the multiplication. + public virtual Matrix Multiply(T scalar) { - if (IsOneT(scalar)) + if (scalar.Equals(One)) { - return; + return Clone(); } - CommonParallel.For( - 0, - RowCount, - i => - { - for (var j = 0; j < ColumnCount; j++) - { - At(i, j, MultiplyT(At(i, j), scalar)); - } - }); + if (scalar.Equals(0.0)) + { + return CreateMatrix(RowCount, ColumnCount); + } + + var result = CreateMatrix(RowCount, ColumnCount); + Multiply(scalar, result); + return result; } /// /// 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 multiply. + /// The matrix to store the result of the multiplication. /// If the result matrix is . /// If the result matrix's dimensions are not the same as this matrix. public virtual void Multiply(T scalar, Matrix result) @@ -146,10 +208,16 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "result"); } - CopyTo(result); - result.Multiply(scalar); + DoMultiply(scalar, result); } + /// + /// 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 abstract void DoMultiply(T scalar, Matrix result); + /// /// Multiplies this matrix by a vector and returns the result. /// @@ -165,7 +233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } /// - /// Multiplies this matrix with a vector and places the results into the result vactor. + /// Multiplies this matrix with a vector and places the results into the result vector. /// /// The vector to multiply with. /// The result of the multiplication. @@ -203,22 +271,17 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CommonParallel.For( - 0, - RowCount, - i => - { - var s = default(T); - for (var j = 0; j != ColumnCount; j++) - { - s = AddT(s, MultiplyT(At(i, j), rightSide[j])); - } - - result[i] = s; - }); + DoMultiply(rightSide, result); } } + /// + /// 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 abstract void DoMultiply(Vector rightSide, Vector result); + /// /// Left multiply a matrix with a vector ( = vector * matrix ). /// @@ -272,22 +335,17 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CommonParallel.For( - 0, - RowCount, - j => - { - var s = default(T); - for (var i = 0; i != leftSide.Count; i++) - { - s = AddT(s, MultiplyT(leftSide[i], At(i, j))); - } - - result[j] = s; - }); + DoLeftMultiply(leftSide, result); } } + /// + /// Left multiply a matrix with a vector ( = vector * matrix ) and place the result in the result vector. + /// + /// The vector to multiply with. + /// The result of the multiplication. + protected abstract void DoLeftMultiply(Vector leftSide, Vector result); + /// /// Multiplies this matrix with another matrix and places the results into the result matrix. /// @@ -327,22 +385,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CommonParallel.For( - 0, - RowCount, - j => - { - for (var i = 0; i != other.ColumnCount; i++) - { - var s = default(T); - for (var l = 0; l < ColumnCount; l++) - { - s = AddT(s, MultiplyT(At(j, l), other.At(l, i))); - } - - result.At(j, i, s); - } - }); + DoMultiply(other, result); } } @@ -370,6 +413,13 @@ namespace MathNet.Numerics.LinearAlgebra.Generic return 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 abstract void DoMultiply(Matrix other, Matrix result); + /// /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix. /// @@ -409,22 +459,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CommonParallel.For( - 0, - RowCount, - j => - { - for (var i = 0; i < RowCount; i++) - { - var s = default(T); - for (var l = 0; l < ColumnCount; l++) - { - s = AddT(s, MultiplyT(At(i, l), other.At(j, l))); - } - - result.At(i, j, AddT(s, result.At(i, j))); - } - }); + DoTransposeAndMultiply(other, result); } } @@ -452,12 +487,23 @@ namespace MathNet.Numerics.LinearAlgebra.Generic return 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 abstract void DoTransposeAndMultiply(Matrix other, Matrix result); + /// /// Negate each element of this matrix. /// - /// If the result matrix is . - /// if the result matrix's dimensions are not the same as this matrix. - public abstract void Negate(); + /// A matrix containing the negated values. + public virtual Matrix Negate() + { + var result = CreateMatrix(RowCount, ColumnCount); + Negate(result); + return result; + } /// /// Negate each element of this matrix and place the results into the result matrix. @@ -477,10 +523,15 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentMatrixDimensions); } - CopyTo(result); - result.Negate(); + DoNegate(result); } + /// + /// Negate each element of this matrix and place the results into the result matrix. + /// + /// The result of the negation. + protected abstract void DoNegate(Matrix result); + /// /// Adds two matrices together and returns the results. /// @@ -739,16 +790,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); } - CommonParallel.For( - 0, - ColumnCount, - j => - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, MultiplyT(At(i, j), other.At(i, j))); - } - }); + DoPointwiseMultiply(other, result); } /// @@ -775,6 +817,13 @@ namespace MathNet.Numerics.LinearAlgebra.Generic return result; } + /// + /// 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 abstract void DoPointwiseMultiply(Matrix other, Matrix result); + /// /// Pointwise divide this matrix by another matrix and stores the result into the result matrix. /// @@ -806,18 +855,16 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); } - CommonParallel.For( - 0, - ColumnCount, - j => - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, DivideT(At(i, j), other.At(i, j))); - } - }); + DoPointwiseDivide(other, result); } + /// + /// 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 abstract void DoPointwiseDivide(Matrix other, Matrix result); + /// /// Generates matrix with random elements. /// @@ -849,17 +896,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// /// The trace of this matrix /// If the matrix is not square - public virtual T Trace() - { - if (RowCount != ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSquare); - } - - var sum = default(T); - CommonParallel.For(0, RowCount, i => sum = AddT(sum, this[i, i])); - return sum; - } + public abstract T Trace(); /// /// Calculates the rank of the matrix diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.cs index 0bd87fe3..036d6411 100644 --- a/src/Numerics/LinearAlgebra/Generic/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Generic/Matrix.cs @@ -38,7 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// /// Defines the base class for Matrix classes. /// - /// Supported data types are double, single, , and . + /// Supported data types are double, single, , and . [Serializable] public abstract partial class Matrix : #if SILVERLIGHT @@ -119,7 +119,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// /// The column of the element. /// - /// The double value to get or set. + /// The value to get or set. /// This method is ranged checked. and /// to get and set values without range checking. public virtual T this[int row, int column] @@ -1492,25 +1492,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// Returns the conjugate transpose of this matrix. /// /// The conjugate transpose of this matrix. - public virtual Matrix ConjugateTranspose() - { - // In case of real return regulart transpose - if (typeof(T) == typeof(double) || (typeof(T) == typeof(float))) - { - return Transpose(); - } - - var ret = CreateMatrix(ColumnCount, RowCount); - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - ret.At(j, i, ConjugateT(At(i, j))); - } - } - - return ret; - } + public abstract Matrix ConjugateTranspose(); /// /// Permute the rows of a matrix according to a permutation. @@ -1788,177 +1770,23 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// Calculates the L1 norm. /// The L1 norm of the matrix. - public virtual double L1Norm() - { - double norm = 0.0; - for (var j = 0; j < ColumnCount; j++) - { - var s = 0.0; - for (var i = 0; i < RowCount; i++) - { - s += AbsoluteT(At(i, j)); - } - - norm = Math.Max(norm, s); - } - - return norm; - } + public abstract T L1Norm(); /// Calculates the L2 norm. /// The L2 norm of the matrix. /// For sparse matrices, the L2 norm is computed using a dense implementation of singular value decomposition. /// In a later release, it will be replaced with a sparse implementation. - public virtual double L2Norm() + public virtual T L2Norm() { return Svd.Create(this, false).Norm2; } /// Calculates the Frobenius norm of this matrix. /// The Frobenius norm of this matrix. - public virtual double FrobeniusNorm() - { - var transpose = Transpose(); - var aat = this * transpose; - - var norm = 0.0; - for (var i = 0; i < RowCount; i++) - { - norm += AbsoluteT(aat.At(i, i)); - } - - norm = Math.Sqrt(norm); - - return norm; - } + public abstract T FrobeniusNorm(); /// Calculates the infinity norm of this matrix. /// The infinity norm of this matrix. - public virtual double InfinityNorm() - { - var norm = 0.0; - for (var i = 0; i < RowCount; i++) - { - var s = 0.0; - for (var j = 0; j < ColumnCount; j++) - { - s += AbsoluteT(At(i, j)); - } - - norm = Math.Max(norm, s); - } - - return norm; - } - - #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); - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected abstract T SubtractT(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); - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected abstract T DivideT(T val1, T val2); - - /// - /// Take absolute value - /// - /// Source value - /// True if one; otherwise false - protected abstract double AbsoluteT(T val1); - - /// - /// Is equal to one? - /// - /// Value to check - /// True if one; otherwise false - private static bool IsOneT(T val1) - { - if (typeof(T) == typeof(Complex)) - { - object obj1 = val1; - return Complex.One.AlmostEqual((Complex)obj1); - } - - if (typeof(T) == typeof(Complex32)) - { - object obj1 = val1; - return Complex32.One.AlmostEqual((Complex32)obj1); - } - - if (typeof(T) == typeof(double)) - { - object obj1 = val1; - return 1.0.AlmostEqualInDecimalPlaces((double)obj1, 15); - } - - if (typeof(T) == typeof(float)) - { - object obj1 = val1; - return 1.0f.AlmostEqualInDecimalPlaces((float)obj1, 7); - } - - throw new NotSupportedException(); - } - - /// - /// Conjugate complex value. In real case the same value is returned - /// - /// Value to conjugate - /// Conjugated value (complex) or the same (real) - private static T ConjugateT(T val1) - { - if (typeof(T) == typeof(Complex)) - { - object obj = val1; - object conj = Complex.Conjugate((Complex)obj); - return (T)conj; - } - - if (typeof(T) == typeof(Complex32)) - { - object obj = val1; - object conj = ((Complex32)obj).Conjugate(); - return (T)conj; - } - - if (typeof(T) == typeof(double)) - { - return val1; - } - - if (typeof(T) == typeof(float)) - { - return val1; - } - - throw new NotSupportedException(); - } - #endregion + public abstract T InfinityNorm(); } } diff --git a/src/Numerics/LinearAlgebra/Generic/Vector.cs b/src/Numerics/LinearAlgebra/Generic/Vector.cs index ef8fb22b..bf4ffc4a 100644 --- a/src/Numerics/LinearAlgebra/Generic/Vector.cs +++ b/src/Numerics/LinearAlgebra/Generic/Vector.cs @@ -55,9 +55,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic private static readonly T Zero = default(T); /// - /// The value on 1.0 for type T. + /// The value of 1.0 for type T. /// - private static readonly T One = SetOne(); + private static readonly T One = Common.SetOne(); /// /// Initializes a new instance of the Vector class. @@ -143,7 +143,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - Add(scalar, result); + DoAdd(scalar, result); return result; } @@ -174,11 +174,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } - DoAdd(scalar, result); } @@ -233,7 +228,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - Add(other, result); + DoAdd(other, result); return result; } @@ -270,15 +265,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = Add(other); - tmp.CopyTo(result); - } - else - { - DoAdd(other, result); - } + DoAdd(other, result); } /// @@ -307,7 +294,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - Subtract(scalar, result); + DoSubtract(scalar, result); return result; } @@ -338,11 +325,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } - DoSubtract(scalar, result); } @@ -394,7 +376,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - Subtract(other, result); + DoSubtract(other, result); return result; } @@ -431,15 +413,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = Subtract(other); - tmp.CopyTo(result); - } - else - { - DoSubtract(other, result); - } + DoSubtract(other, result); } /// @@ -468,7 +442,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - Multiply(scalar, result); + DoMultiply(scalar, result); return result; } @@ -499,11 +473,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } - DoMultiply(scalar, result); } @@ -574,7 +543,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - Divide(scalar, result); + DoDivide(scalar, result); return result; } @@ -605,11 +574,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } - DoDivide(scalar, result); } @@ -644,7 +608,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - PointwiseMultiply(other, result); + DoPointwiseMultiply(other, result); return result; } @@ -679,15 +643,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = PointwiseMultiply(other); - tmp.CopyTo(result); - } - else - { - DoPointwiseMultiply(other, result); - } + DoPointwiseMultiply(other, result); } /// @@ -717,7 +673,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } var result = CreateVector(Count); - PointwiseDivide(other, result); + DoPointwiseDivide(other, result); return result; } @@ -752,15 +708,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = PointwiseDivide(other); - tmp.CopyTo(result); - } - else - { - DoPointwiseDivide(other, result); - } + DoPointwiseDivide(other, result); } /// @@ -829,7 +777,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// /// Matrix M[i,j] = this[i] * v[j]. /// - /// + /// public Matrix OuterProduct(Vector v) { return OuterProduct(this, v); @@ -1586,34 +1534,5 @@ namespace MathNet.Numerics.LinearAlgebra.Generic { CommonParallel.For(0, Count, index => this[index] = default(T)); } - - /// - /// Sets the value of 1.0 for type T. - /// - /// The value of 1.0 for type T. - private static T SetOne() - { - if (typeof(T) == typeof(Complex)) - { - return (T)(object)Complex.One; - } - - if (typeof(T) == typeof(Complex32)) - { - return (T)(object)Complex32.One; - } - - if (typeof(T) == typeof(double)) - { - return (T)(object)1.0; - } - - if (typeof(T) == typeof(float)) - { - return (T)(object)1.0f; - } - - throw new NotSupportedException(); - } } } diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 148d38b1..fd86edac 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -127,19 +127,22 @@ + + Code + + + - - @@ -168,11 +171,7 @@ - - - - @@ -201,11 +200,8 @@ - - - @@ -216,9 +212,7 @@ - - @@ -249,7 +243,6 @@ - @@ -304,7 +297,6 @@ -