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 @@
-