diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
new file mode 100644
index 00000000..5367bf76
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -0,0 +1,742 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// Copyright (c) 2009-2010 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+ using System;
+ using System.Numerics;
+ using Distributions;
+ using Generic;
+ using Properties;
+ using Threading;
+
+ ///
+ /// A Matrix class with dense storage. The underlying storage is a one dimensional array in column-major order.
+ ///
+ public class DenseMatrix : Matrix
+ {
+ ///
+ /// Initializes a new instance of the class. This matrix is square with a given size.
+ ///
+ /// the size of the square matrix.
+ ///
+ /// If is less than one.
+ ///
+ public DenseMatrix(int order)
+ : base(order)
+ {
+ Data = new Complex[order * order];
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ /// The number of rows.
+ ///
+ ///
+ /// The number of columns.
+ ///
+ public DenseMatrix(int rows, int columns)
+ : base(rows, columns)
+ {
+ Data = new Complex[rows * columns];
+ }
+
+ ///
+ /// Initializes a new instance of the class with all entries set to a particular value.
+ ///
+ ///
+ /// The number of rows.
+ ///
+ ///
+ /// The number of columns.
+ ///
+ /// The value which we assign to each element of the matrix.
+ public DenseMatrix(int rows, int columns, Complex value)
+ : base(rows, columns)
+ {
+ Data = new Complex[rows * columns];
+ 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 number of rows.
+ /// The number of columns.
+ /// The one dimensional array to create this matrix from. This array should store the matrix in column-major order.
+ public DenseMatrix(int rows, int columns, Complex[] array)
+ : base(rows, columns)
+ {
+ Data = array;
+ }
+
+ ///
+ /// Initializes a new instance of the class from a 2D array. This constructor
+ /// will allocate a completely new memory block for storing the dense matrix.
+ ///
+ /// The 2D array to create this matrix from.
+ public DenseMatrix(Complex[,] array)
+ : base(array.GetLength(0), array.GetLength(1))
+ {
+ var rows = array.GetLength(0);
+ var columns = array.GetLength(1);
+ Data = new Complex[rows * columns];
+ for (var i = 0; i < rows; i++)
+ {
+ for (var j = 0; j < columns; j++)
+ {
+ Data[(j * rows) + i] = array[i, j];
+ }
+ }
+ }
+
+ ///
+ /// Gets the matrix's data.
+ ///
+ /// The matrix's data.
+ internal Complex[] Data
+ {
+ get;
+ private set;
+ }
+
+ ///
+ /// Creates a DenseMatrix for the given number of rows and columns.
+ ///
+ ///
+ /// The number of rows.
+ ///
+ ///
+ /// The number of columns.
+ ///
+ ///
+ /// A DenseMatrix with the given dimensions.
+ ///
+ public override Matrix CreateMatrix(int numberOfRows, int numberOfColumns)
+ {
+ return new DenseMatrix(numberOfRows, numberOfColumns);
+ }
+
+ ///
+ /// Creates a with a the given dimension.
+ ///
+ /// The size of the vector.
+ ///
+ /// A with the given dimension.
+ ///
+ public override Vector CreateVector(int size)
+ {
+ return new DenseVector(size);
+ }
+
+ ///
+ /// Retrieves the requested element without range checking.
+ ///
+ ///
+ /// The row of the element.
+ ///
+ ///
+ /// The column of the element.
+ ///
+ ///
+ /// The requested element.
+ ///
+ public override Complex At(int row, int column)
+ {
+ return Data[(column * RowCount) + row];
+ }
+
+ ///
+ /// Sets the value of the given element.
+ ///
+ ///
+ /// The row of the element.
+ ///
+ ///
+ /// The column of the element.
+ ///
+ ///
+ /// The value to set the element to.
+ ///
+ public override void At(int row, int column, Complex value)
+ {
+ Data[(column * RowCount) + row] = value;
+ }
+
+ ///
+ /// Sets all values to zero.
+ ///
+ public override void Clear()
+ {
+ Array.Clear(Data, 0, Data.Length);
+ }
+
+ ///
+ /// Returns the transpose of this matrix.
+ ///
+ /// The transpose of this matrix.
+ public override Matrix Transpose()
+ {
+ var ret = new DenseMatrix(ColumnCount, RowCount);
+ for (var j = 0; j < ColumnCount; j++)
+ {
+ var index = j * RowCount;
+ for (var i = 0; i < RowCount; i++)
+ {
+ ret.Data[(i * ColumnCount) + j] = Data[index + i];
+ }
+ }
+
+ return ret;
+ }
+
+ ///
+ /// Returns the conjugate transpose of this matrix.
+ ///
+ /// The conjugate transpose of this matrix.
+ public override Matrix ConjugateTranspose()
+ {
+ var ret = new DenseMatrix(ColumnCount, RowCount);
+ for (var j = 0; j < ColumnCount; j++)
+ {
+ var index = j * RowCount;
+ for (var i = 0; i < RowCount; i++)
+ {
+ ret.Data[(i * ColumnCount) + j] = Data[index + i].Conjugate();
+ }
+ }
+
+ return ret;
+ }
+
+ /// Calculates the L1 norm.
+ /// The L1 norm of the matrix.
+ public override double L1Norm()
+ {
+ var norm = 0.0;
+ for (var j = 0; j < ColumnCount; j++)
+ {
+ var s = 0.0;
+ for (var i = 0; i < RowCount; i++)
+ {
+ s += Data[(j * RowCount) + i].Magnitude;
+ }
+
+ norm = Math.Max(norm, s);
+ }
+
+ return norm;
+ }
+
+ /// Calculates the Frobenius norm of this matrix.
+ /// The Frobenius norm of this matrix.
+ public override double FrobeniusNorm()
+ {
+ var transpose = (DenseMatrix)Transpose();
+ var aat = this * transpose;
+
+ var norm = 0.0;
+ for (var i = 0; i < RowCount; i++)
+ {
+ norm += aat.Data[(i * RowCount) + i].Magnitude;
+ }
+
+ norm = Math.Sqrt(norm);
+ return norm;
+ }
+
+ /// Calculates the infinity norm of this matrix.
+ /// The infinity norm of this matrix.
+ public override 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 += Data[(j * RowCount) + i].Magnitude;
+ }
+
+ norm = Math.Max(norm, s);
+ }
+
+ return norm;
+ }
+
+ #region Elementary operations
+
+ ///
+ /// Adds another matrix to this matrix. The result will be written into this matrix.
+ ///
+ /// The matrix to add to this matrix.
+ /// If the other matrix is .
+ /// If the two matrices don't have the same dimensions.
+ public override void Add(Matrix other)
+ {
+ var m = other as DenseMatrix;
+ if (m == null)
+ {
+ base.Add(other);
+ }
+ else
+ {
+ Add(m);
+ }
+ }
+
+ ///
+ /// 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.
+ ///
+ /// 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 complex.
+ ///
+ /// The complex to multiply with.
+ public override void Multiply(Complex complex)
+ {
+ Control.LinearAlgebraProvider.ScaleArray(complex, 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)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (result == 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);
+ }
+ 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);
+ }
+
+ 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.
+
+ ///
+ /// Initializes a square with all zero's except for ones on the diagonal.
+ ///
+ /// the size of the square matrix.
+ /// A dense identity matrix.
+ ///
+ /// If is less than one.
+ ///
+ public static DenseMatrix Identity(int order)
+ {
+ var m = new DenseMatrix(order);
+ for (var i = 0; i < order; i++)
+ {
+ m[i, i] = Complex.One;
+ }
+
+ return m;
+ }
+
+ #endregion
+
+ ///
+ /// 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 override void Negate()
+ {
+ Multiply(-1);
+ }
+
+ ///
+ /// Generates matrix with random elements.
+ ///
+ /// Number of rows.
+ /// Number of columns.
+ /// Continuous Random Distribution or Source
+ ///
+ /// An numberOfRows-by-numberOfColumns matrix with elements distributed according to the provided distribution.
+ ///
+ /// If the parameter is not positive.
+ /// If the parameter is not positive.
+ public override Matrix Random(int numberOfRows, int numberOfColumns, IContinuousDistribution distribution)
+ {
+ if (numberOfRows < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfRows");
+ }
+
+ if (numberOfColumns < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfColumns");
+ }
+
+ var matrix = CreateMatrix(numberOfRows, numberOfColumns);
+ CommonParallel.For(
+ 0,
+ ColumnCount,
+ j =>
+ {
+ for (var i = 0; i < matrix.RowCount; i++)
+ {
+ matrix[i, j] = new Complex(distribution.Sample(), distribution.Sample());
+ }
+ });
+
+ return matrix;
+ }
+
+ ///
+ /// Generates matrix with random elements.
+ ///
+ /// Number of rows.
+ /// Number of columns.
+ /// Continuous Random Distribution or Source
+ ///
+ /// An numberOfRows-by-numberOfColumns matrix with elements distributed according to the provided distribution.
+ ///
+ /// If the parameter is not positive.
+ /// If the parameter is not positive.
+ public override Matrix Random(int numberOfRows, int numberOfColumns, IDiscreteDistribution distribution)
+ {
+ if (numberOfRows < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfRows");
+ }
+
+ if (numberOfColumns < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfColumns");
+ }
+
+ var matrix = CreateMatrix(numberOfRows, numberOfColumns);
+ CommonParallel.For(
+ 0,
+ ColumnCount,
+ j =>
+ {
+ for (var i = 0; i < matrix.RowCount; i++)
+ {
+ matrix[i, j] = new Complex(distribution.Sample(), distribution.Sample());
+ }
+ });
+
+ return matrix;
+ }
+
+ #region Simple arithmetic of type T
+ ///
+ /// Add two values T+T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of addition
+ protected sealed override Complex AddT(Complex val1, Complex val2)
+ {
+ return val1 + val2;
+ }
+
+ ///
+ /// Subtract two values T-T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of subtract
+ protected sealed override Complex SubtractT(Complex val1, Complex val2)
+ {
+ return val1 - val2;
+ }
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Divide two values T/T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of divide
+ protected sealed override Complex DivideT(Complex val1, Complex val2)
+ {
+ return val1 / val2;
+ }
+
+ ///
+ /// Is equal to one?
+ ///
+ /// Value to check
+ /// True if one; otherwise false
+ protected sealed override bool IsOneT(Complex val1)
+ {
+ return val1.AlmostEqual(Complex.One);
+ }
+
+ ///
+ /// Take absolute value
+ ///
+ /// Source alue
+ /// True if one; otherwise false
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs
new file mode 100644
index 00000000..52f9c5ba
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs
@@ -0,0 +1,1513 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// Copyright (c) 2009-2010 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+ using System;
+ using System.Collections.Generic;
+ using System.Numerics;
+ using Distributions;
+ using Generic;
+ using NumberTheory;
+ using Properties;
+ using Threading;
+
+ ///
+ /// A vector using dense storage.
+ ///
+ public class DenseVector : Vector
+ {
+ ///
+ /// Initializes a new instance of the class with a given size.
+ ///
+ ///
+ /// the size of the vector.
+ ///
+ ///
+ /// If is less than one.
+ ///
+ public DenseVector(int size) : base(size)
+ {
+ Data = new Complex[size];
+ }
+
+ ///
+ /// Initializes a new instance of the class with a given size
+ /// and each element set to the given value;
+ ///
+ ///
+ /// the size of the vector.
+ ///
+ ///
+ /// the value to set each element to.
+ ///
+ ///
+ /// If is less than one.
+ ///
+ public DenseVector(int size, Complex value) : this(size)
+ {
+ for (var index = 0; index < Data.Length; index++)
+ {
+ Data[index] = value;
+ }
+ }
+
+ ///
+ /// Initializes a new instance of the class by
+ /// copying the values from another.
+ ///
+ ///
+ /// The vector to create the new vector from.
+ ///
+ public DenseVector(Vector other) : this(other.Count)
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => this[index] = other[index]);
+ }
+
+ ///
+ /// Initializes a new instance of the class by
+ /// copying the values from another.
+ ///
+ ///
+ /// The vector to create the new vector from.
+ ///
+ public DenseVector(DenseVector other) : this(other.Count)
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => Data[index] = other.Data[index]);
+ }
+
+ ///
+ /// Initializes a new instance of the class for an array.
+ ///
+ /// The array to create this vector from.
+ /// The vector does not copy the array, but keeps a reference to it. Any
+ /// changes to the vector will also change the array.
+ public DenseVector(Complex[] array) : base(array.Length)
+ {
+ Data = array;
+ }
+
+ ///
+ /// Gets the vector's internal data.
+ ///
+ /// The vector's internal data.
+ /// Changing values in the array also changes the corresponding value in vector. Use with care.
+ internal Complex[] Data
+ {
+ get;
+ private set;
+ }
+
+ ///
+ /// Returns a reference to the internal data structure.
+ ///
+ /// The DenseVector whose internal data we are
+ /// returning.
+ ///
+ /// A reference to the internal date of the given vector.
+ ///
+ public static implicit operator Complex[](DenseVector vector)
+ {
+ if (vector == null)
+ {
+ throw new ArgumentNullException();
+ }
+
+ return vector.Data;
+ }
+
+ ///
+ /// Returns a vector bound directly to a reference of the provided array.
+ ///
+ /// The array to bind to the DenseVector object.
+ ///
+ /// A DenseVector whose values are bound to the given array.
+ ///
+ public static implicit operator DenseVector(Complex[] array)
+ {
+ if (array == null)
+ {
+ throw new ArgumentNullException();
+ }
+
+ return new DenseVector(array);
+ }
+
+ ///
+ /// Create a matrix based on this vector in column form (one single column).
+ ///
+ /// This vector as a column matrix.
+ public override Matrix ToColumnMatrix()
+ {
+ var matrix = new DenseMatrix(Count, 1);
+ for (var i = 0; i < Data.Length; i++)
+ {
+ matrix[i, 0] = Data[i];
+ }
+
+ return matrix;
+ }
+
+ ///
+ /// Create a matrix based on this vector in row form (one single row).
+ ///
+ /// This vector as a row matrix.
+ public override Matrix ToRowMatrix()
+ {
+ var matrix = new DenseMatrix(1, Count);
+ for (var i = 0; i < Data.Length; i++)
+ {
+ matrix[0, i] = Data[i];
+ }
+
+ return matrix;
+ }
+
+ /// Gets or sets the value at the given .
+ /// The index of the value to get or set.
+ /// The value of the vector at the given .
+ /// If is negative or
+ /// greater than the size of the vector.
+ public override Complex this[int index]
+ {
+ get
+ {
+ return Data[index];
+ }
+
+ set
+ {
+ Data[index] = value;
+ }
+ }
+
+ ///
+ /// Creates a matrix with the given dimensions using the same storage type
+ /// as this vector.
+ ///
+ ///
+ /// The number of rows.
+ ///
+ ///
+ /// The number of columns.
+ ///
+ ///
+ /// A matrix with the given dimensions.
+ ///
+ public override Matrix CreateMatrix(int rows, int columns)
+ {
+ return new DenseMatrix(rows, columns);
+ }
+
+ ///
+ /// Creates a Vector of the given size using the same storage type
+ /// as this vector.
+ ///
+ ///
+ /// The size of the Vector to create.
+ ///
+ ///
+ /// The new Vector.
+ ///
+ public override Vector CreateVector(int size)
+ {
+ return new DenseVector(size);
+ }
+
+ ///
+ /// Adds a complex to each element of the vector.
+ ///
+ /// The complex to add.
+ /// A copy of the vector with the complex added.
+ public override Vector Add(Complex complex)
+ {
+ if (complex == Complex.Zero)
+ {
+ return Clone();
+ }
+
+ var copy = (DenseVector)Clone();
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => copy.Data[index] += complex);
+ return copy;
+ }
+
+ ///
+ /// Adds a complex to each element of the vector and stores the result in the result vector.
+ ///
+ /// The complex to add.
+ /// The vector to store the result of the addition.
+ /// If the result vector is .
+ /// If this vector and are not the same size.
+ public override void Add(Complex complex, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ var dense = result as DenseVector;
+ if (dense == null)
+ {
+ base.Add(complex, result);
+ }
+ else
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => dense.Data[index] = Data[index] + complex);
+ }
+ }
+
+ ///
+ /// Adds another vector to this vector.
+ ///
+ /// The vector to add to this one.
+ /// A new vector containing the sum of both vectors.
+ /// If the other vector is .
+ /// If this vector and are not the same size.
+ public override Vector Add(Vector other)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ var denseVector = other as DenseVector;
+
+ if (denseVector == null)
+ {
+ return base.Add(other);
+ }
+
+ var copy = (DenseVector)Clone();
+ Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, Complex.One, denseVector.Data);
+ return copy;
+ }
+
+ ///
+ /// Adds another vector to this vector and stores the result into the result vector.
+ ///
+ /// The vector to add to this one.
+ /// The vector to store the result of the addition.
+ /// If the other vector is .
+ /// If the result vector is .
+ /// If this vector and are not the same size.
+ /// If this vector and are not the same size.
+ public override void Add(Vector other, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ if (Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
+ {
+ var tmp = Add(other);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ var rdense = result as DenseVector;
+ var odense = other as DenseVector;
+ if (rdense != null && odense != null)
+ {
+ CopyTo(result);
+ Control.LinearAlgebraProvider.AddVectorToScaledVector(rdense.Data, Complex.One, odense.Data);
+ }
+ else
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => result[index] = Data[index] + other[index]);
+ }
+ }
+ }
+
+ ///
+ /// Returns a Vector containing the same values of .
+ ///
+ /// This method is included for completeness.
+ /// The vector to get the values from.
+ /// A vector containing a the same values as .
+ /// If is .
+ public static Vector operator +(DenseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException("rightSide");
+ }
+
+ return rightSide.Plus();
+ }
+
+ ///
+ /// Adds two Vectors together and returns the results.
+ ///
+ /// One of the vectors to add.
+ /// The other vector to add.
+ /// The result of the addition.
+ /// If and are not the same size.
+ /// If or is .
+ public static Vector operator +(DenseVector leftSide, DenseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException("rightSide");
+ }
+
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException("leftSide");
+ }
+
+ if (leftSide.Count != rightSide.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rightSide");
+ }
+
+ return leftSide.Add(rightSide);
+ }
+
+ ///
+ /// Subtracts a complex from each element of the vector.
+ ///
+ /// The complex to subtract.
+ /// A new vector containing the subtraction of this vector and the complex.
+ public override Vector Subtract(Complex complex)
+ {
+ if (complex == Complex.Zero)
+ {
+ return Clone();
+ }
+
+ var copy = (DenseVector)Clone();
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => copy.Data[index] -= complex);
+ return copy;
+ }
+
+ ///
+ /// Subtracts a complex from each element of the vector and stores the result in the result vector.
+ ///
+ /// The complex to subtract.
+ /// The vector to store the result of the subtraction.
+ /// If the result vector is .
+ /// If this vector and are not the same size.
+ public override void Subtract(Complex complex, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ var dense = result as DenseVector;
+ if (dense == null)
+ {
+ base.Subtract(complex, result);
+ }
+ else
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => dense.Data[index] = Data[index] - complex);
+ }
+ }
+
+ ///
+ /// Subtracts another vector from this vector.
+ ///
+ /// The vector to subtract from this one.
+ /// A new vector containing the subtraction of the the two vectors.
+ /// If the other vector is .
+ /// If this vector and are not the same size.
+ public override Vector Subtract(Vector other)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ var denseVector = other as DenseVector;
+
+ if (denseVector == null)
+ {
+ return base.Subtract(other);
+ }
+
+ var copy = (DenseVector)Clone();
+ Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, -Complex.One, denseVector.Data);
+ return copy;
+ }
+
+ ///
+ /// Subtracts another vector to this vector and stores the result into the result vector.
+ ///
+ /// The vector to subtract from this one.
+ /// The vector to store the result of the subtraction.
+ /// If the other vector is .
+ /// If the result vector is .
+ /// If this vector and are not the same size.
+ /// If this vector and are not the same size.
+ public override void Subtract(Vector other, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ if (Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
+ {
+ var tmp = Subtract(other);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ var rdense = result as DenseVector;
+ var odense = other as DenseVector;
+ if (rdense != null && odense != null)
+ {
+ CopyTo(result);
+ Control.LinearAlgebraProvider.AddVectorToScaledVector(rdense.Data, -Complex.One, odense.Data);
+ }
+ else
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => result[index] = Data[index] - other[index]);
+ }
+ }
+ }
+
+ ///
+ /// Returns a Vector containing the negated values of .
+ ///
+ /// The vector to get the values from.
+ /// A vector containing the negated values as .
+ /// If is .
+ public static Vector operator -(DenseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException("rightSide");
+ }
+
+ return rightSide.Negate();
+ }
+
+ ///
+ /// Subtracts two Vectors and returns the results.
+ ///
+ /// The vector to subtract from.
+ /// The vector to subtract.
+ /// The result of the subtraction.
+ /// If and are not the same size.
+ /// If or is .
+ public static Vector operator -(DenseVector leftSide, DenseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException("rightSide");
+ }
+
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException("leftSide");
+ }
+
+ if (leftSide.Count != rightSide.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rightSide");
+ }
+
+ return leftSide.Subtract(rightSide);
+ }
+
+ ///
+ /// Returns a negated vector.
+ ///
+ /// The negated vector.
+ /// Added as an alternative to the unary negation operator.
+ public override Vector Negate()
+ {
+ var result = new DenseVector(Count);
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => result[index] = -Data[index]);
+
+ return result;
+ }
+
+ ///
+ /// Multiplies a complex to each element of the vector.
+ ///
+ /// The complex to multiply.
+ /// A new vector that is the multiplication of the vector and the complex.
+ public override Vector Multiply(Complex complex)
+ {
+ if (complex == Complex.One)
+ {
+ return Clone();
+ }
+
+ var copy = (DenseVector)Clone();
+ Control.LinearAlgebraProvider.ScaleArray(complex, copy.Data);
+ return copy;
+ }
+
+ ///
+ /// Computes the dot product between this vector and another vector.
+ ///
+ /// The other vector to add.
+ /// The result of the addition.
+ /// If is not of the same size.
+ /// If is .
+ public override Complex DotProduct(Vector other)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ var denseVector = other as DenseVector;
+
+ if (denseVector == null)
+ {
+ return base.DotProduct(other);
+ }
+
+ return Control.LinearAlgebraProvider.DotProduct(Data, denseVector.Data);
+ }
+
+ ///
+ /// Multiplies a vector with a complex.
+ ///
+ /// The vector to scale.
+ /// The Complex value.
+ /// The result of the multiplication.
+ /// If is .
+ public static DenseVector operator *(DenseVector leftSide, Complex rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException("leftSide");
+ }
+
+ return (DenseVector)leftSide.Multiply(rightSide);
+ }
+
+ ///
+ /// Multiplies a vector with a complex.
+ ///
+ /// The Complex value.
+ /// The vector to scale.
+ /// The result of the multiplication.
+ /// If is .
+ public static DenseVector operator *(Complex leftSide, DenseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException("rightSide");
+ }
+
+ return (DenseVector)rightSide.Multiply(leftSide);
+ }
+
+ ///
+ /// Computes the dot product between two Vectors.
+ ///
+ /// The left row vector.
+ /// The right column vector.
+ /// The dot product between the two vectors.
+ /// If and are not the same size.
+ /// If or is .
+ public static Complex operator *(DenseVector leftSide, DenseVector rightSide)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException("rightSide");
+ }
+
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException("leftSide");
+ }
+
+ if (leftSide.Count != rightSide.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rightSide");
+ }
+
+ return Control.LinearAlgebraProvider.DotProduct(leftSide.Data, rightSide.Data);
+ }
+
+ ///
+ /// Divides a vector with a complex.
+ ///
+ /// The vector to divide.
+ /// The Complex value.
+ /// The result of the division.
+ /// If is .
+ public static DenseVector operator /(DenseVector leftSide, Complex rightSide)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException("leftSide");
+ }
+
+ return (DenseVector)leftSide.Multiply(Complex.One / rightSide);
+ }
+
+ ///
+ /// Returns the index of the absolute minimum element.
+ ///
+ /// The index of absolute minimum element.
+ public override int AbsoluteMinimumIndex()
+ {
+ var index = 0;
+ var min = Data[index].Magnitude;
+ for (var i = 1; i < Count; i++)
+ {
+ var test = Data[i].Magnitude;
+ if (test < min)
+ {
+ index = i;
+ min = test;
+ }
+ }
+
+ return index;
+ }
+
+ ///
+ /// Returns the value of the absolute minimum element.
+ ///
+ /// The value of the absolute minimum element.
+ public override double AbsoluteMinimum()
+ {
+ return Data[AbsoluteMinimumIndex()].Magnitude;
+ }
+
+ ///
+ /// Returns the value of the absolute maximum element.
+ ///
+ /// The value of the absolute maximum element.
+ public override double AbsoluteMaximum()
+ {
+ return Data[AbsoluteMaximumIndex()].Magnitude;
+ }
+
+ ///
+ /// Returns the index of the absolute maximum element.
+ ///
+ /// The index of absolute maximum element.
+ public override int AbsoluteMaximumIndex()
+ {
+ var index = 0;
+ var max = Data[index].Magnitude;
+ for (var i = 1; i < Count; i++)
+ {
+ var test = Data[i].Magnitude;
+ if (test > max)
+ {
+ index = i;
+ max = test;
+ }
+ }
+
+ return index;
+ }
+
+ ///
+ /// Creates a vector containing specified elements.
+ ///
+ /// The first element to begin copying from.
+ /// The number of elements to copy.
+ /// A vector containing a copy of the specified elements.
+ /// - If is not positive or
+ /// greater than or equal to the size of the vector.
+ /// - If + is greater than or equal to the size of the vector.
+ ///
+ /// If is not positive.
+ public override Vector SubVector(int index, int length)
+ {
+ if (index < 0 || index >= Count)
+ {
+ throw new ArgumentOutOfRangeException("index");
+ }
+
+ if (length <= 0)
+ {
+ throw new ArgumentOutOfRangeException("length");
+ }
+
+ if (index + length > Count)
+ {
+ throw new ArgumentOutOfRangeException("length");
+ }
+
+ var result = new DenseVector(length);
+
+ CommonParallel.For(
+ index,
+ index + length,
+ i => result.Data[i - index] = Data[i]);
+ return result;
+ }
+
+ ///
+ /// Set the values of this vector to the given values.
+ ///
+ /// The array containing the values to use.
+ /// If is .
+ /// If is not the same size as this vector.
+ public override void SetValues(Complex[] values)
+ {
+ if (values == null)
+ {
+ throw new ArgumentNullException("values");
+ }
+
+ if (values.Length != Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "values");
+ }
+
+ CommonParallel.For(
+ 0,
+ values.Length,
+ i => Data[i] = values[i]);
+ }
+
+ ///
+ /// Computes the sum of the vector's elements.
+ ///
+ /// The sum of the vector's elements.
+ public override Complex Sum()
+ {
+ var result = Complex.Zero;
+ for (var i = 0; i < Count; i++)
+ {
+ result += Data[i];
+ }
+
+ return result;
+ }
+
+ ///
+ /// Computes the sum of the absolute value of the vector's elements.
+ ///
+ /// The sum of the absolute value of the vector's elements.
+ public override double SumMagnitudes()
+ {
+ double result = 0;
+ for (var i = 0; i < Count; i++)
+ {
+ result += Data[i].Magnitude;
+ }
+
+ return result;
+ }
+
+ ///
+ /// Pointwise multiplies this vector with another vector.
+ ///
+ /// The vector to pointwise multiply with this one.
+ /// A new vector which is the pointwise multiplication of the two vectors.
+ /// If the other vector is .
+ /// If this vector and are not the same size.
+ public override Vector PointwiseMultiply(Vector other)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ var denseVector = other as DenseVector;
+
+ if (denseVector == null)
+ {
+ return base.PointwiseMultiply(other);
+ }
+
+ var copy = (DenseVector)Clone();
+ CommonParallel.For(
+ 0,
+ Count,
+ index => copy[index] *= other[index]);
+ return copy;
+ }
+
+ ///
+ /// Pointwise multiplies this vector with another vector and stores the result into the result vector.
+ ///
+ /// The vector to pointwise multiply with this one.
+ /// The vector to store the result of the pointwise multiplication.
+ /// If the other vector is .
+ /// If the result vector is .
+ /// If this vector and are not the same size.
+ /// If this vector and are not the same size.
+ public override void PointwiseMultiply(Vector other, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ if (Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
+ {
+ var tmp = PointwiseMultiply(other);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ var dense = result as DenseVector;
+ if (dense == null)
+ {
+ base.PointwiseMultiply(other, result);
+ }
+ else
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => dense.Data[index] = Data[index] * other[index]);
+ }
+ }
+ }
+
+ ///
+ /// Pointwise divide this vector with another vector.
+ ///
+ /// The vector to pointwise divide this one by.
+ /// A new vector which is the pointwise division of the two vectors.
+ /// If the other vector is .
+ /// If this vector and are not the same size.
+ public override Vector PointwiseDivide(Vector other)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ var denseVector = other as DenseVector;
+
+ if (denseVector == null)
+ {
+ return base.PointwiseMultiply(other);
+ }
+
+ var copy = (DenseVector)Clone();
+ CommonParallel.For(
+ 0,
+ Count,
+ index => copy[index] /= other[index]);
+ return copy;
+ }
+
+ ///
+ /// Pointwise divide this vector with another vector and stores the result into the result vector.
+ ///
+ /// The vector to pointwise divide this one by.
+ /// The vector to store the result of the pointwise division.
+ /// If the other vector is .
+ /// If the result vector is .
+ /// If this vector and are not the same size.
+ /// If this vector and are not the same size.
+ public override void PointwiseDivide(Vector other, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (Count != other.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+ }
+
+ if (Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
+ {
+ var tmp = PointwiseDivide(other);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ var dense = result as DenseVector;
+ if (dense == null)
+ {
+ base.PointwiseDivide(other, result);
+ }
+ else
+ {
+ CommonParallel.For(
+ 0,
+ Data.Length,
+ index => dense.Data[index] = Data[index] / other[index]);
+ }
+ }
+ }
+
+ ///
+ /// Outer product of two vectors
+ ///
+ /// First vector
+ /// Second vector
+ /// Matrix M[i,j] = u[i]*v[j]
+ /// If the u vector is .
+ /// If the v vector is .
+ public static DenseMatrix OuterProduct(DenseVector u, DenseVector v)
+ {
+ if (u == null)
+ {
+ throw new ArgumentNullException("u");
+ }
+
+ if (v == null)
+ {
+ throw new ArgumentNullException("v");
+ }
+
+ var matrix = new DenseMatrix(u.Count, v.Count);
+ CommonParallel.For(
+ 0,
+ u.Count,
+ i =>
+ {
+ for (var j = 0; j < v.Count; j++)
+ {
+ matrix.At(i, j, u.Data[i] * v.Data[j]);
+ }
+ });
+ return matrix;
+ }
+
+ ///
+ /// Generates a vector with random elements
+ ///
+ /// Number of elements in the vector.
+ /// Continuous Random Distribution or Source
+ ///
+ /// A vector with n-random elements distributed according
+ /// to the specified random distribution.
+ ///
+ /// If the n vector is non positive.
+ public override Vector Random(int length, IContinuousDistribution randomDistribution)
+ {
+ if (length < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
+ }
+
+ var v = (DenseVector)CreateVector(length);
+ for (var index = 0; index < v.Data.Length; index++)
+ {
+ v.Data[index] = randomDistribution.Sample();
+ }
+
+ return v;
+ }
+
+ ///
+ /// Generates a vector with random elements
+ ///
+ /// Number of elements in the vector.
+ /// Continuous Random Distribution or Source
+ ///
+ /// A vector with n-random elements distributed according
+ /// to the specified random distribution.
+ ///
+ /// If the n vector is non positive.
+ public override Vector Random(int length, IDiscreteDistribution randomDistribution)
+ {
+ if (length < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
+ }
+
+ var v = (DenseVector)CreateVector(length);
+ for (var index = 0; index < v.Data.Length; index++)
+ {
+ v.Data[index] = randomDistribution.Sample();
+ }
+
+ return v;
+ }
+
+ ///
+ /// Tensor Product (Dyadic) of this and another vector.
+ ///
+ /// The vector to operate on.
+ ///
+ /// Matrix M[i,j] = this[i] * v[j].
+ ///
+ ///
+ public Matrix TensorMultiply(DenseVector v)
+ {
+ return OuterProduct(this, v);
+ }
+
+ ///
+ /// Computes the p-Norm.
+ ///
+ /// The p value.
+ /// Scalar ret = (sum(abs(this[i])^p))^(1/p)
+ public override double Norm(double p)
+ {
+ if (p < 0.0)
+ {
+ throw new ArgumentOutOfRangeException("p");
+ }
+
+ if (1.0 == p)
+ {
+ return CommonParallel.Aggregate(
+ 0,
+ Count,
+ index => Data[index].Magnitude);
+ }
+
+ if (Double.IsPositiveInfinity(p))
+ {
+ return CommonParallel.Select(
+ 0,
+ Count,
+ (index, localData) => Math.Max(localData, Data[index].Magnitude),
+ Math.Max);
+ }
+
+ var sum = CommonParallel.Aggregate(
+ 0,
+ Count,
+ index => Math.Pow(Data[index].Magnitude, p));
+
+ return Math.Pow(sum, 1.0 / p);
+ }
+
+ ///
+ /// Normalizes this vector to a unit vector with respect to the p-norm.
+ ///
+ ///
+ /// The p value.
+ ///
+ ///
+ /// This vector normalized to a unit vector with respect to the p-norm.
+ ///
+ public override Vector Normalize(double p)
+ {
+ if (p < 0.0)
+ {
+ throw new ArgumentOutOfRangeException("p");
+ }
+
+ var norm = Norm(p);
+ var clone = Clone();
+ if (norm == 0.0)
+ {
+ return clone;
+ }
+
+ clone.Multiply(1.0 / norm, clone);
+
+ return clone;
+ }
+
+ #region Parse Functions
+
+ ///
+ /// Creates a Complex dense vector based on a string. The string can be in the following formats (without the
+ /// quotes): 'n', 'n;n;..', '(n;n;..)', '[n;n;...]', where n is a Complex.
+ ///
+ ///
+ /// A Complex dense vector containing the values specified by the given string.
+ ///
+ ///
+ /// The string to parse.
+ ///
+ public static DenseVector Parse(string value)
+ {
+ return Parse(value, null);
+ }
+
+ ///
+ /// Creates a Complex dense vector based on a string. The string can be in the following formats (without the
+ /// quotes): 'n', 'n;n;..', '(n;n;..)', '[n;n;...]', where n is a double.
+ ///
+ ///
+ /// A Complex dense vector containing the values specified by the given string.
+ ///
+ ///
+ /// the string to parse.
+ ///
+ ///
+ /// An that supplies culture-specific formatting information.
+ ///
+ public static DenseVector Parse(string value, IFormatProvider formatProvider)
+ {
+ if (value == null)
+ {
+ throw new ArgumentNullException(value);
+ }
+
+ value = value.Trim();
+ if (value.Length == 0)
+ {
+ throw new FormatException();
+ }
+
+ // strip out parens
+ if (value.StartsWith("(", StringComparison.Ordinal))
+ {
+ if (!value.EndsWith(")", StringComparison.Ordinal))
+ {
+ throw new FormatException();
+ }
+
+ value = value.Substring(1, value.Length - 2).Trim();
+ }
+
+ if (value.StartsWith("[", StringComparison.Ordinal))
+ {
+ if (!value.EndsWith("]", StringComparison.Ordinal))
+ {
+ throw new FormatException();
+ }
+
+ value = value.Substring(1, value.Length - 2).Trim();
+ }
+
+ // keywords
+ var textInfo = formatProvider.GetTextInfo();
+ var keywords = new[] { textInfo.ListSeparator };
+
+ // lexing
+ var tokens = new LinkedList();
+ GlobalizationHelper.Tokenize(tokens.AddFirst(value), keywords, 0);
+ var token = tokens.First;
+
+ if (token == null || tokens.Count.IsEven())
+ {
+ throw new FormatException();
+ }
+
+ // parsing
+ var data = new Complex[(tokens.Count + 1) >> 1];
+ for (var i = 0; i < data.Length; i++)
+ {
+ if (token == null || token.Value == textInfo.ListSeparator)
+ {
+ throw new FormatException();
+ }
+
+ data[i] = token.Value.ToComplex(formatProvider);
+
+ token = token.Next;
+ if (token != null)
+ {
+ token = token.Next;
+ }
+ }
+
+ return new DenseVector(data);
+ }
+
+ ///
+ /// Converts the string representation of a complex dense vector to double-precision dense vector equivalent.
+ /// A return value indicates whether the conversion succeeded or failed.
+ ///
+ ///
+ /// A string containing a complex vector to convert.
+ ///
+ ///
+ /// The parsed value.
+ ///
+ ///
+ /// If the conversion succeeds, the result will contain a complex number equivalent to value.
+ /// Otherwise the result will be null.
+ ///
+ public static bool TryParse(string value, out DenseVector result)
+ {
+ return TryParse(value, null, out result);
+ }
+
+ ///
+ /// Converts the string representation of a complex dense vector to double-precision dense vector equivalent.
+ /// A return value indicates whether the conversion succeeded or failed.
+ ///
+ ///
+ /// A string containing a complex vector to convert.
+ ///
+ ///
+ /// An that supplies culture-specific formatting information about value.
+ ///
+ ///
+ /// The parsed value.
+ ///
+ ///
+ /// If the conversion succeeds, the result will contain a complex number equivalent to value.
+ /// Otherwise the result will be null.
+ ///
+ public static bool TryParse(string value, IFormatProvider formatProvider, out DenseVector result)
+ {
+ bool ret;
+ try
+ {
+ result = Parse(value, formatProvider);
+ ret = true;
+ }
+ catch (ArgumentNullException)
+ {
+ result = null;
+ ret = false;
+ }
+ catch (FormatException)
+ {
+ result = null;
+ ret = false;
+ }
+
+ return ret;
+ }
+
+ #endregion
+
+ ///
+ /// Returns the index of the absolute maximum element.
+ ///
+ /// The index of absolute maximum element.
+ public override int MaximumIndex()
+ {
+ throw new NotSupportedException();
+ }
+
+ ///
+ /// Returns the index of the minimum element.
+ ///
+ /// The index of minimum element.
+ public override int MinimumIndex()
+ {
+ throw new NotSupportedException();
+ }
+
+ ///
+ /// Resets all values to zero.
+ ///
+ public override void Clear()
+ {
+ Array.Clear(Data, 0, Data.Length);
+ }
+
+ ///
+ /// Conjugates vector and save result to
+ ///
+ /// Target vector
+ public override void Conjugate(Vector target)
+ {
+ if (target == null)
+ {
+ throw new ArgumentNullException("target");
+ }
+
+ if (Count != target.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "target");
+ }
+
+ if (ReferenceEquals(this, target))
+ {
+ var tmp = CreateVector(Count);
+ Conjugate(tmp);
+ tmp.CopyTo(target);
+ }
+
+ CommonParallel.For(
+ 0,
+ Count,
+ index => target[index] = this[index].Conjugate());
+ }
+
+ #region Simple arithmetic of type T
+ ///
+ /// Add two values T+T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of addition
+ protected sealed override Complex AddT(Complex val1, Complex val2)
+ {
+ return val1 + val2;
+ }
+
+ ///
+ /// Subtract two values T-T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of subtract
+ protected sealed override Complex SubtractT(Complex val1, Complex val2)
+ {
+ return val1 - val2;
+ }
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Divide two values T/T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of divide
+ protected sealed override Complex DivideT(Complex val1, Complex val2)
+ {
+ return val1 / val2;
+ }
+
+ ///
+ /// Is equal to one?
+ ///
+ /// Value to check
+ /// True if one; otherwise false
+ protected sealed override bool IsOneT(Complex val1)
+ {
+ return Complex.One.AlmostEqual(val1);
+ }
+
+ ///
+ /// Take absolute value
+ ///
+ /// Source alue
+ /// True if one; otherwise false
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
new file mode 100644
index 00000000..0bf51c79
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -0,0 +1,1718 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// Copyright (c) 2009-2010 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+ using System;
+ using System.Linq;
+ using System.Numerics;
+ using Distributions;
+ using Generic;
+ using Properties;
+ using Threading;
+
+ ///
+ /// A matrix type for diagonal matrices.
+ ///
+ ///
+ /// Diagonal matrices can be non-square matrices but the diagonal always starts
+ /// at element 0,0. A diagonal matrix will throw an exception if non diagonal
+ /// entries are set. The exception to this is when the off diagonal elements are
+ /// 0.0 or NaN; these settings will cause no change to the diagonal matrix.
+ ///
+ public class DiagonalMatrix : Matrix
+ {
+ ///
+ /// Initializes a new instance of the class. This matrix is square with a given size.
+ ///
+ /// the size of the square matrix.
+ ///
+ /// If is less than one.
+ ///
+ public DiagonalMatrix(int order) : base(order)
+ {
+ Data = new Complex[order * order];
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ /// The number of rows.
+ ///
+ ///
+ /// The number of columns.
+ ///
+ public DiagonalMatrix(int rows, int columns) : base(rows, columns)
+ {
+ Data = new Complex[Math.Min(rows, columns)];
+ }
+
+ ///
+ /// Initializes a new instance of the class with all entries set to a particular value.
+ ///
+ ///
+ /// The number of rows.
+ ///
+ ///
+ /// The number of columns.
+ ///
+ /// The value which we assign to each element of the matrix.
+ public DiagonalMatrix(int rows, int columns, Complex value) : base(rows, columns)
+ {
+ Data = new Complex[Math.Min(rows, columns)];
+ for (var i = 0; i < Data.Length; i++)
+ {
+ Data[i] = value;
+ }
+ }
+
+ ///
+ /// Initializes a new instance of the class from a one dimensional array with diagonal elements. This constructor
+ /// will reference the one dimensional array and not copy it.
+ ///
+ /// The number of rows.
+ /// The number of columns.
+ /// The one dimensional array which contain diagonal elements.
+ public DiagonalMatrix(int rows, int columns, Complex[] diagonalArray) : base(rows, columns)
+ {
+ Data = diagonalArray;
+ }
+
+ ///
+ /// Initializes a new instance of the class from a 2D array.
+ ///
+ /// The 2D array to create this matrix from.
+ /// When contains an off-diagonal element.
+ /// Depending on the implementation, an
+ /// may be thrown if one of the indices is outside the dimensions of the matrix.
+ public DiagonalMatrix(Complex[,] array) : this(array.GetLength(0), array.GetLength(1))
+ {
+ var rows = array.GetLength(0);
+ var columns = array.GetLength(1);
+
+ for (var i = 0; i < rows; i++)
+ {
+ for (var j = 0; j < columns; j++)
+ {
+ if (i == j)
+ {
+ Data[i] = array[i, j];
+ }
+ else if (((array[i, j].Real != 0.0) && !double.IsNaN(array[i, j].Real)) || ((array[i, j].Imaginary != 0.0) && !double.IsNaN(array[i, j].Imaginary)))
+ {
+ throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
+ }
+ }
+ }
+ }
+
+ ///
+ /// Gets the matrix's data.
+ ///
+ /// The matrix's data.
+ internal Complex[] Data
+ {
+ get;
+ private set;
+ }
+
+ ///
+ /// Retrieves the requested element without range checking.
+ ///
+ ///
+ /// The row of the element.
+ ///
+ ///
+ /// The column of the element.
+ ///
+ ///
+ /// The requested element.
+ ///
+ /// Depending on the implementation, an
+ /// may be thrown if one of the indices is outside the dimensions of the matrix.
+ public override Complex At(int row, int column)
+ {
+ return row == column ? Data[row] : 0.0;
+ }
+
+ ///
+ /// Sets the value of the given element.
+ ///
+ ///
+ /// The row of the element.
+ ///
+ ///
+ /// The column of the element.
+ ///
+ ///
+ /// The value to set the element to.
+ ///
+ /// When trying to set an off diagonal element.
+ /// Depending on the implementation, an
+ /// may be thrown if one of the indices is outside the dimensions of the matrix.
+ public override void At(int row, int column, Complex value)
+ {
+ if (row == column)
+ {
+ Data[row] = value;
+ }
+ else if (((value.Real != 0.0) && !double.IsNaN(value.Real)) || ((value.Imaginary != 0.0) && !double.IsNaN(value.Imaginary)))
+ {
+ throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
+ }
+ }
+
+ ///
+ /// Creates a DiagonalMatrix for the given number of rows and columns.
+ ///
+ ///
+ /// The number of rows.
+ ///
+ ///
+ /// The number of columns.
+ ///
+ ///
+ /// A DiagonalMatrix with the given dimensions.
+ ///
+ public override Matrix CreateMatrix(int numberOfRows, int numberOfColumns)
+ {
+ return new DiagonalMatrix(numberOfRows, numberOfColumns);
+ }
+
+ ///
+ /// Creates a with a the given dimension.
+ ///
+ /// The size of the vector.
+ ///
+ /// A with the given dimension.
+ ///
+ public override Vector CreateVector(int size)
+ {
+ return new SparseVector(size);
+ }
+
+ ///
+ /// Sets all values to zero.
+ ///
+ public override void Clear()
+ {
+ Array.Clear(Data, 0, Data.Length);
+ }
+
+ ///
+ /// Indicates whether the current object is equal to another object of the same type.
+ ///
+ ///
+ /// An object to compare with this object.
+ ///
+ ///
+ /// true if the current object is equal to the parameter; otherwise, false.
+ ///
+ public override bool Equals(object obj)
+ {
+ var diagonalMatrix = obj as DiagonalMatrix;
+
+ if (diagonalMatrix == null)
+ {
+ return base.Equals(obj);
+ }
+
+ // Accept if the argument is the same object as this
+ if (ReferenceEquals(this, diagonalMatrix))
+ {
+ return true;
+ }
+
+ if (diagonalMatrix.Data.Length != Data.Length)
+ {
+ return false;
+ }
+
+ // If all else fails, perform element wise comparison.
+ return !Data.Where((t, i) => t != diagonalMatrix.Data[i]).Any();
+ }
+
+ ///
+ /// Returns a hash code for this instance.
+ ///
+ ///
+ /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table.
+ ///
+ public override int GetHashCode()
+ {
+ var hashNum = Math.Min(Data.Length, 25);
+ long hash = 0;
+ for (var i = 0; i < hashNum; i++)
+ {
+#if SILVERLIGHT
+ hash ^= Precision.DoubleToInt64Bits(Data[i].GetHashCode());
+#else
+ hash ^= BitConverter.DoubleToInt64Bits(Data[i].GetHashCode());
+#endif
+ }
+
+ return BitConverter.ToInt32(BitConverter.GetBytes(hash), 4);
+ }
+
+ #region Elementary operations
+ ///
+ /// Adds another matrix to this matrix. The result will be written into this matrix.
+ ///
+ /// The matrix to add to this matrix.
+ /// If the other matrix is .
+ /// If the two matrices don't have the same dimensions.
+ /// If is not .
+ public override void Add(Matrix other)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ var m = other as DiagonalMatrix;
+ if (m == null)
+ {
+ throw new ArgumentException(Resources.ArgumentTypeMismatch);
+ }
+
+ Add(m);
+ }
+
+ ///
+ /// 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(DiagonalMatrix 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.
+ ///
+ /// The matrix to subtract.
+ /// If the other matrix is .
+ /// If the two matrices don't have the same dimensions.
+ /// If is not .
+ public override void Subtract(Matrix other)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ var m = other as DiagonalMatrix;
+ if (m == null)
+ {
+ throw new ArgumentException(Resources.ArgumentTypeMismatch);
+ }
+
+ 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(DiagonalMatrix 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);
+ }
+
+ ///
+ /// Copies the values of the given array to the diagonal.
+ ///
+ /// The array to copy the values from. The length of the vector should be
+ /// Min(Rows, Columns).
+ /// If is .
+ /// If the length of does not
+ /// equal Min(Rows, Columns).
+ /// For non-square matrices, the elements of are copied to
+ /// this[i,i].
+ public override void SetDiagonal(Complex[] source)
+ {
+ if (source == null)
+ {
+ throw new ArgumentNullException("source");
+ }
+
+ if (source.Length != Data.Length)
+ {
+ throw new ArgumentException(Resources.ArgumentArraysSameLength, "source");
+ }
+
+ CommonParallel.For(0, source.Length, index => Data[index] = source[index]);
+ }
+
+ ///
+ /// Copies the values of the given to the diagonal.
+ ///
+ /// The vector to copy the values from. The length of the vector should be
+ /// Min(Rows, Columns).
+ /// If is .
+ /// If the length of does not
+ /// equal Min(Rows, Columns).
+ /// For non-square matrices, the elements of are copied to
+ /// this[i,i].
+ public override void SetDiagonal(Vector source)
+ {
+ var denseSource = source as DenseVector;
+ if (denseSource == null)
+ {
+ base.SetDiagonal(source);
+ return;
+ }
+
+ if (Data.Length != denseSource.Data.Length)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "source");
+ }
+
+ CommonParallel.For(0, denseSource.Data.Length, index => Data[index] = denseSource.Data[index]);
+ }
+
+ ///
+ /// Multiplies each element of this matrix with a scalar.
+ ///
+ /// The scalar to multiply with.
+ public override void Multiply(Complex scalar)
+ {
+ if (scalar == 0.0)
+ {
+ Clear();
+ return;
+ }
+
+ if (scalar == 1.0)
+ {
+ return;
+ }
+
+ Control.LinearAlgebraProvider.ScaleArray(scalar, Data);
+ }
+
+ ///
+ /// Multiplies this diagonal matrix with another diagonal matrix and places the results into the result diagonal 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)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (result == 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 DiagonalMatrix;
+ var r = result as DiagonalMatrix;
+
+ if (m == null || r == null)
+ {
+ base.Multiply(other, result);
+ }
+ else
+ {
+ var thisDataCopy = new Complex[r.Data.Length];
+ var otherDataCopy = new Complex[r.Data.Length];
+
+ CommonParallel.For(0, (r.Data.Length > Data.Length) ? Data.Length : r.Data.Length, index => thisDataCopy[index] = Data[index]);
+ CommonParallel.For(0, (r.Data.Length > m.Data.Length) ? m.Data.Length : r.Data.Length, index => otherDataCopy[index] = m.Data[index]);
+
+ Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, 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 DiagonalMatrix;
+ if (m == null)
+ {
+ return base.Multiply(other);
+ }
+
+ var result = (DiagonalMatrix)CreateMatrix(RowCount, other.ColumnCount);
+ Multiply(other, result);
+ return result;
+ }
+
+ ///
+ /// Multiplies this matrix with a vector and places the results into the result matrix.
+ ///
+ /// The vector to multiply with.
+ /// The result of the multiplication.
+ /// If is .
+ /// If is .
+ /// If result.Count != this.RowCount.
+ /// If this.ColumnCount != .Count.
+ public override void Multiply(Vector rightSide, Vector result)
+ {
+ if (rightSide == null)
+ {
+ throw new ArgumentNullException("rightSide");
+ }
+
+ if (ColumnCount != rightSide.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "rightSide");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (RowCount != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ if (ReferenceEquals(rightSide, result))
+ {
+ var tmp = result.CreateVector(result.Count);
+ Multiply(rightSide, tmp);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ // Clear the result vector
+ result.Clear();
+
+ // Multiply the elements in the vector with the corresponding diagonal element in this.
+ for (var r = 0; r < Data.Length; r++)
+ {
+ result[r] = Data[r] * rightSide[r];
+ }
+ }
+ }
+
+ ///
+ /// 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.
+ /// If is .
+ /// If the result matrix is .
+ /// If result.Count != this.ColumnCount.
+ /// If this.RowCount != .Count.
+ public override void LeftMultiply(Vector leftSide, Vector result)
+ {
+ if (leftSide == null)
+ {
+ throw new ArgumentNullException("leftSide");
+ }
+
+ if (RowCount != leftSide.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "leftSide");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (ColumnCount != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ if (ReferenceEquals(leftSide, result))
+ {
+ var tmp = result.CreateVector(result.Count);
+ LeftMultiply(leftSide, tmp);
+ tmp.CopyTo(result);
+ }
+ else
+ {
+ // Clear the result vector
+ result.Clear();
+
+ // Multiply the elements in the vector with the corresponding diagonal element in this.
+ for (var r = 0; r < Data.Length; r++)
+ {
+ result[r] = Data[r] * leftSide[r];
+ }
+ }
+ }
+
+ ///
+ /// Computes the determinant of this matrix.
+ ///
+ /// The determinant of this matrix.
+ public override Complex Determinant()
+ {
+ if (RowCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare);
+ }
+
+ return Data.Aggregate(Complex.One, (current, t) => current * t);
+ }
+
+ ///
+ /// Returns the elements of the diagonal in a .
+ ///
+ /// The elements of the diagonal.
+ /// For non-square matrices, the method returns Min(Rows, Columns) elements where
+ /// i == j (i is the row index, and j is the column index).
+ public override Vector Diagonal()
+ {
+ // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
+ // not influence onto diagonal elements
+ return new DenseVector((Complex[])Data.Clone());
+ }
+
+ ///
+ /// Multiplies this diagonal matrix with transpose of another diagonal matrix and places the results into the result diagonal 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 otherDiagonal = other as DiagonalMatrix;
+ var resultDiagonal = result as DiagonalMatrix;
+
+ if (otherDiagonal == null || resultDiagonal == null)
+ {
+ base.TransposeAndMultiply(other, result);
+ return;
+ }
+
+ Multiply(otherDiagonal.Transpose(), result);
+ }
+
+ ///
+ /// 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 otherDiagonal = other as DiagonalMatrix;
+ if (otherDiagonal == null)
+ {
+ return base.TransposeAndMultiply(other);
+ }
+
+ if (ColumnCount != otherDiagonal.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var result = (DiagonalMatrix)CreateMatrix(RowCount, other.RowCount);
+ TransposeAndMultiply(other, result);
+ return result;
+ }
+
+ ///
+ /// Multiplies two diagonal 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 DiagonalMatrix operator *(DiagonalMatrix leftSide, DiagonalMatrix 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 (DiagonalMatrix)leftSide.Multiply(rightSide);
+ }
+
+ #endregion
+
+ ///
+ /// Copies the elements of this matrix to the given matrix.
+ ///
+ ///
+ /// The matrix to copy values into.
+ ///
+ ///
+ /// If target is .
+ ///
+ ///
+ /// If this and the target matrix do not have the same dimensions..
+ ///
+ public override void CopyTo(Matrix target)
+ {
+ var diagonalTarget = target as DiagonalMatrix;
+
+ if (diagonalTarget == null)
+ {
+ base.CopyTo(target);
+ return;
+ }
+
+ if (ReferenceEquals(this, target))
+ {
+ return;
+ }
+
+ if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "target");
+ }
+
+ CommonParallel.For(0, Data.Length, index => diagonalTarget.Data[index] = Data[index]);
+ }
+
+ ///
+ /// Returns the transpose of this matrix.
+ ///
+ /// The transpose of this matrix.
+ public override Matrix Transpose()
+ {
+ var ret = new DiagonalMatrix(ColumnCount, RowCount);
+ CommonParallel.For(0, Data.Length, index => ret.Data[index] = Data[index]);
+ return ret;
+ }
+
+ ///
+ /// Returns the conjugate transpose of this matrix.
+ ///
+ /// The conjugate transpose of this matrix.
+ public override Matrix ConjugateTranspose()
+ {
+ var ret = new DiagonalMatrix(ColumnCount, RowCount);
+ CommonParallel.For(0, Data.Length, index => ret.Data[index] = Data[index].Conjugate());
+ return ret;
+ }
+
+ ///
+ /// Copies the requested column elements into the given vector.
+ ///
+ /// The column to copy elements from.
+ /// The row to start copying from.
+ /// The number of elements to copy.
+ /// The to copy the column into.
+ /// If the result is .
+ /// If is negative,
+ /// or greater than or equal to the number of columns.
+ /// If is negative,
+ /// or greater than or equal to the number of rows.
+ /// If +
+ /// is greater than or equal to the number of rows.
+ /// If is not positive.
+ /// If result.Count < length.
+ public override void Column(int columnIndex, int rowIndex, int length, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (columnIndex >= ColumnCount || columnIndex < 0)
+ {
+ throw new ArgumentOutOfRangeException("columnIndex");
+ }
+
+ if (rowIndex >= RowCount || rowIndex < 0)
+ {
+ throw new ArgumentOutOfRangeException("rowIndex");
+ }
+
+ if (rowIndex + length > RowCount)
+ {
+ throw new ArgumentOutOfRangeException("length");
+ }
+
+ if (length < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
+ }
+
+ if (result.Count < length)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ // Clear the result and copy the diagonal entry.
+ result.Clear();
+ if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < Data.Length)
+ {
+ result[columnIndex - rowIndex] = Data[columnIndex];
+ }
+ }
+
+ ///
+ /// Copies the requested row elements into a new .
+ ///
+ /// The row to copy elements from.
+ /// The column to start copying from.
+ /// The number of elements to copy.
+ /// The to copy the column into.
+ /// If the result is .
+ /// If is negative,
+ /// or greater than or equal to the number of columns.
+ /// If is negative,
+ /// or greater than or equal to the number of rows.
+ /// If +
+ /// is greater than or equal to the number of rows.
+ /// If is not positive.
+ /// If result.Count < length.
+ public override void Row(int rowIndex, int columnIndex, int length, Vector result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (rowIndex >= RowCount || rowIndex < 0)
+ {
+ throw new ArgumentOutOfRangeException("rowIndex");
+ }
+
+ if (columnIndex >= ColumnCount || columnIndex < 0)
+ {
+ throw new ArgumentOutOfRangeException("columnIndex");
+ }
+
+ if (columnIndex + length > ColumnCount)
+ {
+ throw new ArgumentOutOfRangeException("length");
+ }
+
+ if (length < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
+ }
+
+ if (result.Count < length)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
+ }
+
+ // Clear the result and copy the diagonal entry.
+ result.Clear();
+ if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < Data.Length)
+ {
+ result[rowIndex - columnIndex] = Data[rowIndex];
+ }
+ }
+
+ /// Calculates the L1 norm.
+ /// The L1 norm of the matrix.
+ public override double L1Norm()
+ {
+ return Data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
+ }
+
+ /// Calculates the L2 norm.
+ /// The L2 norm of the matrix.
+ public override double L2Norm()
+ {
+ return Data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
+ }
+
+ /// Calculates the Frobenius norm of this matrix.
+ /// The Frobenius norm of this matrix.
+ public override double FrobeniusNorm()
+ {
+ var norm = Data.Sum(t => t.Magnitude * t.Magnitude);
+ return Math.Sqrt(norm);
+ }
+
+ /// Calculates the infinity norm of this matrix.
+ /// The infinity norm of this matrix.
+ public override double InfinityNorm()
+ {
+ return L1Norm();
+ }
+
+ /// Calculates the condition number of this matrix.
+ /// The condition number of the matrix.
+ public override double ConditionNumber()
+ {
+ var maxSv = double.NegativeInfinity;
+ var minSv = double.PositiveInfinity;
+ for (var i = 0; i < Data.Length; i++)
+ {
+ maxSv = Math.Max(maxSv, Data[i].Magnitude);
+ minSv = Math.Min(minSv, Data[i].Magnitude);
+ }
+
+ return maxSv / minSv;
+ }
+
+ /// Computes the inverse of this matrix.
+ /// If is not a square matrix.
+ /// If is singular.
+ /// The inverse of this matrix.
+ public override Matrix Inverse()
+ {
+ if (RowCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare);
+ }
+
+ var inverse = (DiagonalMatrix)Clone();
+ for (var i = 0; i < Data.Length; i++)
+ {
+ if (Data[i] != 0.0)
+ {
+ inverse.Data[i] = 1.0 / Data[i];
+ }
+ else
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixNotSingular);
+ }
+ }
+
+ return inverse;
+ }
+
+ ///
+ /// Returns a new matrix containing the lower triangle of this matrix.
+ ///
+ /// The lower triangle of this matrix.
+ public override Matrix LowerTriangle()
+ {
+ return Clone();
+ }
+
+ ///
+ /// Puts the lower triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void LowerTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ if (ReferenceEquals(this, result))
+ {
+ return;
+ }
+
+ result.Clear();
+ for (var i = 0; i < Data.Length; i++)
+ {
+ result[i, i] = Data[i];
+ }
+ }
+
+ ///
+ /// 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()
+ {
+ return new DiagonalMatrix(RowCount, ColumnCount);
+ }
+
+ ///
+ /// Puts the strictly lower triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void StrictlyLowerTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ result.Clear();
+ }
+
+ ///
+ /// Returns a new matrix containing the upper triangle of this matrix.
+ ///
+ /// The upper triangle of this matrix.
+ public override Matrix UpperTriangle()
+ {
+ return Clone();
+ }
+
+ ///
+ /// Puts the upper triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void UpperTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ result.Clear();
+ for (var i = 0; i < Data.Length; i++)
+ {
+ result[i, i] = Data[i];
+ }
+ }
+
+ ///
+ /// 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()
+ {
+ return new DiagonalMatrix(RowCount, ColumnCount);
+ }
+
+ ///
+ /// Puts the strictly upper triangle of this matrix into the result matrix.
+ ///
+ /// Where to store the lower triangle.
+ /// If is .
+ /// If the result matrix's dimensions are not the same as this matrix.
+ public override void StrictlyUpperTriangle(Matrix result)
+ {
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ result.Clear();
+ }
+
+ ///
+ /// Creates a matrix that contains the values from the requested sub-matrix.
+ ///
+ /// The row to start copying from.
+ /// The number of rows to copy. Must be positive.
+ /// The column to start copying from.
+ /// The number of columns to copy. Must be positive.
+ /// The requested sub-matrix.
+ /// If: - is
+ /// negative, or greater than or equal to the number of rows.
+ /// - is negative, or greater than or equal to the number
+ /// of columns.
+ /// - (columnIndex + columnLength) >= Columns
+ /// - (rowIndex + rowLength) >= Rows
+ /// If or
+ /// is not positive.
+ public override Matrix SubMatrix(int rowIndex, int rowLength, int columnIndex, int columnLength)
+ {
+ if (rowIndex >= RowCount || rowIndex < 0)
+ {
+ throw new ArgumentOutOfRangeException("rowIndex");
+ }
+
+ if (columnIndex >= ColumnCount || columnIndex < 0)
+ {
+ throw new ArgumentOutOfRangeException("columnIndex");
+ }
+
+ if (rowLength < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "rowLength");
+ }
+
+ if (columnLength < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "columnLength");
+ }
+
+ var colMax = columnIndex + columnLength;
+ var rowMax = rowIndex + rowLength;
+
+ if (rowMax > RowCount)
+ {
+ throw new ArgumentOutOfRangeException("rowLength");
+ }
+
+ if (colMax > ColumnCount)
+ {
+ throw new ArgumentOutOfRangeException("columnLength");
+ }
+
+ var result = new SparseMatrix(rowLength, columnLength);
+
+ if (rowIndex > columnIndex && columnIndex + columnLength > rowIndex)
+ {
+ for (var i = 0; rowIndex - columnIndex + i < Math.Min(columnLength, rowLength); i++)
+ {
+ result[i, rowIndex - columnIndex + i] = Data[rowIndex + i];
+ }
+ }
+ else if (rowIndex < columnIndex && rowIndex + rowLength > columnIndex)
+ {
+ for (var i = 0; rowIndex - columnIndex + i < Math.Min(columnLength, rowLength); i++)
+ {
+ result[columnIndex - rowIndex + i, i] = Data[columnIndex + i];
+ }
+ }
+ else
+ {
+ for (var i = 0; i < Math.Min(columnLength, rowLength); i++)
+ {
+ result[i, i] = Data[rowIndex + i];
+ }
+ }
+
+ return result;
+ }
+
+ ///
+ /// Returns this matrix as a multidimensional array.
+ ///
+ /// A multidimensional containing the values of this matrix.
+ public override Complex[,] ToArray()
+ {
+ var result = new Complex[RowCount, ColumnCount];
+ for (var i = 0; i < Data.Length; i++)
+ {
+ result[i, i] = Data[i];
+ }
+
+ return result;
+ }
+
+ ///
+ /// Creates a new and inserts the given column at the given index.
+ ///
+ /// The index of where to insert the column.
+ /// The column to insert.
+ /// A new with the inserted column.
+ /// If is .
+ /// If is < zero or > the number of columns.
+ /// If the size of != the number of rows.
+ public override Matrix InsertColumn(int columnIndex, Vector column)
+ {
+ if (column == null)
+ {
+ throw new ArgumentNullException("column");
+ }
+
+ if (columnIndex < 0 || columnIndex > ColumnCount)
+ {
+ throw new ArgumentOutOfRangeException("columnIndex");
+ }
+
+ if (column.Count != RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension, "column");
+ }
+
+ var result = new SparseMatrix(RowCount, ColumnCount + 1);
+
+ for (var i = 0; i < columnIndex; i++)
+ {
+ result.SetColumn(i, Column(i));
+ }
+
+ result.SetColumn(columnIndex, column);
+
+ for (var i = columnIndex + 1; i < ColumnCount + 1; i++)
+ {
+ result.SetColumn(i, Column(i - 1));
+ }
+
+ return result;
+ }
+
+ ///
+ /// Creates a new and inserts the given row at the given index.
+ ///
+ /// The index of where to insert the row.
+ /// The row to insert.
+ /// A new with the inserted column.
+ /// If is .
+ /// If is < zero or > the number of rows.
+ /// If the size of != the number of columns.
+ public override Matrix InsertRow(int rowIndex, Vector row)
+ {
+ if (row == null)
+ {
+ throw new ArgumentNullException("row");
+ }
+
+ if (rowIndex < 0 || rowIndex > RowCount)
+ {
+ throw new ArgumentOutOfRangeException("rowIndex");
+ }
+
+ if (row.Count != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension, "row");
+ }
+
+ var result = new SparseMatrix(RowCount + 1, ColumnCount);
+
+ for (var i = 0; i < rowIndex; i++)
+ {
+ result.SetRow(i, Row(i));
+ }
+
+ result.SetRow(rowIndex, row);
+
+ for (var i = rowIndex + 1; i < RowCount; i++)
+ {
+ result.SetRow(i, Row(i - 1));
+ }
+
+ return result;
+ }
+
+ ///
+ /// Stacks this matrix on top of the given matrix and places the result into the result .
+ ///
+ /// The matrix to stack this matrix upon.
+ /// The combined .
+ /// If lower is .
+ /// If upper.Columns != lower.Columns.
+ public override Matrix Stack(Matrix lower)
+ {
+ if (lower == null)
+ {
+ throw new ArgumentNullException("lower");
+ }
+
+ if (lower.ColumnCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
+ }
+
+ var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount);
+ Stack(lower, result);
+ return result;
+ }
+
+ ///
+ /// Stacks this matrix on top of the given matrix and places the result into the result .
+ ///
+ /// The matrix to stack this matrix upon.
+ /// The combined .
+ /// If lower is .
+ /// If upper.Columns != lower.Columns.
+ public override void Stack(Matrix lower, Matrix result)
+ {
+ if (lower == null)
+ {
+ throw new ArgumentNullException("lower");
+ }
+
+ if (lower.ColumnCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.RowCount != (RowCount + lower.RowCount) || result.ColumnCount != ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ // Clear the result matrix
+ result.Clear();
+
+ // Copy the diagonal part into the result matrix.
+ for (var i = 0; i < Data.Length; i++)
+ {
+ result[i, i] = Data[i];
+ }
+
+ // Copy the lower matrix into the result matrix.
+ for (var i = 0; i < lower.RowCount; i++)
+ {
+ for (var j = 0; j < lower.ColumnCount; j++)
+ {
+ result[i + RowCount, j] = lower[i, j];
+ }
+ }
+ }
+
+ ///
+ /// Concatenates this matrix with the given matrix.
+ ///
+ /// The matrix to concatenate.
+ /// The combined .
+ public override Matrix Append(Matrix right)
+ {
+ if (right == null)
+ {
+ throw new ArgumentNullException("right");
+ }
+
+ if (right.RowCount != RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ var result = new SparseMatrix(RowCount, ColumnCount + right.ColumnCount);
+ Append(right, result);
+ return result;
+ }
+
+ ///
+ /// Concatenates this matrix with the given matrix and places the result into the result .
+ ///
+ /// The matrix to concatenate.
+ /// The combined .
+ public override void Append(Matrix right, Matrix result)
+ {
+ if (right == null)
+ {
+ throw new ArgumentNullException("right");
+ }
+
+ if (right.RowCount != RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.ColumnCount != (ColumnCount + right.ColumnCount) || result.RowCount != RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ // Clear the result matrix
+ result.Clear();
+
+ // Copy the diagonal part into the result matrix.
+ for (var i = 0; i < Data.Length; i++)
+ {
+ result[i, i] = Data[i];
+ }
+
+ // Copy the lower matrix into the result matrix.
+ for (var i = 0; i < right.RowCount; i++)
+ {
+ for (var j = 0; j < right.ColumnCount; j++)
+ {
+ result[i, j + RowCount] = right[i, j];
+ }
+ }
+ }
+
+ ///
+ /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix,
+ /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
+ /// The values of off the off diagonal matrices/blocks are set to zero.
+ ///
+ /// The lower, right matrix.
+ /// If lower is .
+ /// the combined matrix
+ public override Matrix DiagonalStack(Matrix lower)
+ {
+ if (lower == null)
+ {
+ throw new ArgumentNullException("lower");
+ }
+
+ var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
+ DiagonalStack(lower, result);
+ return result;
+ }
+
+ ///
+ /// Diagonally stacks his matrix on top of the given matrix and places the combined matrix into the result matrix.
+ ///
+ /// The lower, right matrix.
+ /// The combined matrix
+ /// If lower is .
+ /// If the result matrix is .
+ /// If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).
+ public override void DiagonalStack(Matrix lower, Matrix result)
+ {
+ if (lower == null)
+ {
+ throw new ArgumentNullException("lower");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.RowCount != RowCount + lower.RowCount || result.ColumnCount != ColumnCount + lower.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ // Clear the result matrix
+ result.Clear();
+
+ // Copy the diagonal part into the result matrix.
+ for (var i = 0; i < Data.Length; i++)
+ {
+ result[i, i] = Data[i];
+ }
+
+ // Copy the lower matrix into the result matrix.
+ CommonParallel.For(0, lower.RowCount, i => CommonParallel.For(0, lower.ColumnCount, j => result.At(i + RowCount, j + ColumnCount, lower.At(i, j))));
+ }
+
+ ///
+ /// 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.
+ /// If the other matrix is .
+ /// If the result matrix is .
+ /// If this matrix and are not the same size.
+ /// If this matrix and are not the same size.
+ public override void PointwiseMultiply(Matrix other, Matrix result)
+ {
+ if (other == null)
+ {
+ throw new ArgumentNullException("other");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (ColumnCount != other.ColumnCount || RowCount != other.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ if (ColumnCount != result.ColumnCount || RowCount != result.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+ }
+
+ var m = other as DiagonalMatrix;
+ var r = result as DiagonalMatrix;
+
+ if (m == null || r == null)
+ {
+ base.PointwiseMultiply(other, result);
+ }
+ else
+ {
+ Control.LinearAlgebraProvider.PointWiseMultiplyArrays(Data, m.Data, r.Data);
+ }
+ }
+
+ ///
+ /// Permute the columns of a matrix according to a permutation.
+ ///
+ /// The column permutation to apply to this matrix.
+ /// Always thrown
+ /// Permutation in diagonal matrix are senseless, because of matrix nature
+ public override void PermuteColumns(Permutation p)
+ {
+ throw new InvalidOperationException("Permutations in diagonal matrix are not allowed");
+ }
+
+ ///
+ /// Permute the rows of a matrix according to a permutation.
+ ///
+ /// The row permutation to apply to this matrix.
+ /// Always thrown
+ /// Permutation in diagonal matrix are senseless, because of matrix nature
+ public override void PermuteRows(Permutation p)
+ {
+ throw new InvalidOperationException("Permutations in diagonal matrix are not allowed");
+ }
+
+ #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 diagonal identity matrix.
+ ///
+ /// If is less than one.
+ ///
+ public static DiagonalMatrix Identity(int order)
+ {
+ var m = new DiagonalMatrix(order);
+ for (var i = 0; i < order; i++)
+ {
+ m.Data[i] = 1.0;
+ }
+
+ return m;
+ }
+
+ #endregion
+
+ ///
+ /// Negates each element of this matrix.
+ ///
+ public override void Negate()
+ {
+ Multiply(-1);
+ }
+
+ ///
+ /// Generates matrix with random elements.
+ ///
+ /// Number of rows.
+ /// Number of columns.
+ /// Continuous Random Distribution or Source
+ ///
+ /// An numberOfRows-by-numberOfColumns matrix with elements distributed according to the provided distribution.
+ ///
+ /// If the parameter is not positive.
+ /// If the parameter is not positive.
+ public override Matrix Random(int numberOfRows, int numberOfColumns, IContinuousDistribution distribution)
+ {
+ if (numberOfRows < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfRows");
+ }
+
+ if (numberOfColumns < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfColumns");
+ }
+
+ var matrix = CreateMatrix(numberOfRows, numberOfColumns);
+ CommonParallel.For(
+ 0,
+ ColumnCount,
+ j =>
+ {
+ for (var i = 0; i < matrix.RowCount; i++)
+ {
+ matrix[i, j] = new Complex(distribution.Sample(), distribution.Sample());
+ }
+ });
+
+ return matrix;
+ }
+
+ ///
+ /// Generates matrix with random elements.
+ ///
+ /// Number of rows.
+ /// Number of columns.
+ /// Continuous Random Distribution or Source
+ ///
+ /// An numberOfRows-by-numberOfColumns matrix with elements distributed according to the provided distribution.
+ ///
+ /// If the parameter is not positive.
+ /// If the parameter is not positive.
+ public override Matrix Random(int numberOfRows, int numberOfColumns, IDiscreteDistribution distribution)
+ {
+ if (numberOfRows < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfRows");
+ }
+
+ if (numberOfColumns < 1)
+ {
+ throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfColumns");
+ }
+
+ var matrix = CreateMatrix(numberOfRows, numberOfColumns);
+ CommonParallel.For(
+ 0,
+ ColumnCount,
+ j =>
+ {
+ for (var i = 0; i < matrix.RowCount; i++)
+ {
+ matrix[i, j] = new Complex(distribution.Sample(), distribution.Sample());
+ }
+ });
+
+ return matrix;
+ }
+
+ #region Simple arithmetic of type T
+ ///
+ /// Add two values T+T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of addition
+ protected sealed override Complex AddT(Complex val1, Complex val2)
+ {
+ return val1 + val2;
+ }
+
+ ///
+ /// Subtract two values T-T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of subtract
+ protected sealed override Complex SubtractT(Complex val1, Complex val2)
+ {
+ return val1 - val2;
+ }
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Divide two values T/T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of divide
+ protected sealed override Complex DivideT(Complex val1, Complex val2)
+ {
+ return val1 / val2;
+ }
+
+ ///
+ /// Is equal to one?
+ ///
+ /// Value to check
+ /// True if one; otherwise false
+ protected sealed override bool IsOneT(Complex val1)
+ {
+ return Complex.One.AlmostEqual(val1);
+ }
+
+ ///
+ /// Take absolute value
+ ///
+ /// Source alue
+ /// True if one; otherwise false
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs
new file mode 100644
index 00000000..374a8a64
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs
@@ -0,0 +1,223 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+ using Threading;
+
+ ///
+ /// A class which encapsulates the functionality of a Cholesky factorization for dense matrices.
+ /// For a symmetric, positive definite matrix A, the Cholesky factorization
+ /// is an lower triangular matrix L so that A = L*L'.
+ ///
+ ///
+ /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
+ /// or positive definite, the constructor will throw an exception.
+ ///
+ public class DenseCholesky : Cholesky
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// Cholesky factorization when the constructor is called and cache it's factorization.
+ ///
+ /// The matrix to factor.
+ /// If is null.
+ /// If is not a square matrix.
+ /// If is not positive definite.
+ public DenseCholesky(DenseMatrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare);
+ }
+
+ // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
+ var factor = (DenseMatrix)matrix.Clone();
+ Control.LinearAlgebraProvider.CholeskyFactor(factor.Data, factor.RowCount);
+ CholeskyFactor = factor;
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A Cholesky factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (result.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ if (result.ColumnCount != input.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ if (input.RowCount != CholeskyFactor.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var dinput = input as DenseMatrix;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment.");
+ }
+
+ var dresult = result as DenseMatrix;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment.");
+ }
+
+ // Copy the contents of input to result.
+ CommonParallel.For(0, dinput.Data.Length, index => dresult.Data[index] = dinput.Data[index]);
+
+ // Cholesky solve by overwriting result.
+ var dfactor = (DenseMatrix)CholeskyFactor;
+ Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount);
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A Cholesky factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (input.Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ if (input.Count != CholeskyFactor.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var dinput = input as DenseVector;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment.");
+ }
+
+ var dresult = result as DenseVector;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment.");
+ }
+
+ // Copy the contents of input to result.
+ CommonParallel.For(0, dinput.Data.Length, index => dresult.Data[index] = dinput.Data[index]);
+
+ // Cholesky solve by overwriting result.
+ var dfactor = (DenseMatrix)CholeskyFactor;
+ Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1);
+ }
+
+ #region Simple arithmetic of type T
+ ///
+ /// Add two values T+T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of addition
+ protected sealed override Complex AddT(Complex val1, Complex val2)
+ {
+ return val1 + val2;
+ }
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Returns the natural (base e) logarithm of a specified number.
+ ///
+ /// A number whose logarithm is to be found
+ /// Natural (base e) logarithm
+ protected sealed override Complex LogT(Complex val1)
+ {
+ return val1.NaturalLogarithm();
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
new file mode 100644
index 00000000..baf323fb
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
@@ -0,0 +1,224 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+ using Threading;
+
+ ///
+ /// A class which encapsulates the functionality of an LU factorization.
+ /// For a matrix A, the LU factorization is a pair of lower triangular matrix L and
+ /// upper triangular matrix U so that A = L*U.
+ ///
+ ///
+ /// The computation of the LU factorization is done at construction time.
+ ///
+ public class DenseLU : LU
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// LU factorization when the constructor is called and cache it's factorization.
+ ///
+ /// The matrix to factor.
+ /// If is null.
+ /// If is not a square matrix.
+ public DenseLU(DenseMatrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare);
+ }
+
+ // Create an array for the pivot indices.
+ Pivots = new int[matrix.RowCount];
+
+ // Create a new matrix for the LU factors, then perform factorization (while overwriting).
+ var factors = (DenseMatrix)matrix.Clone();
+ Control.LinearAlgebraProvider.LUFactor(factors.Data, factors.RowCount, Pivots);
+ Factors = factors;
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A LU factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (result.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ if (result.ColumnCount != input.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ if (input.RowCount != Factors.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var dinput = input as DenseMatrix;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment.");
+ }
+
+ var dresult = result as DenseMatrix;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment.");
+ }
+
+ // Copy the contents of input to result.
+ CommonParallel.For(0, dinput.Data.Length, index => dresult.Data[index] = dinput.Data[index]);
+
+ // LU solve by overwriting result.
+ var dfactors = (DenseMatrix)Factors;
+ Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A LU factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (input.Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ if (input.Count != Factors.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var dinput = input as DenseVector;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment.");
+ }
+
+ var dresult = result as DenseVector;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment.");
+ }
+
+ // Copy the contents of input to result.
+ CommonParallel.For(0, dinput.Data.Length, index => dresult.Data[index] = dinput.Data[index]);
+
+ // LU solve by overwriting result.
+ var dfactors = (DenseMatrix)Factors;
+ Control.LinearAlgebraProvider.LUSolveFactored(1, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data);
+ }
+
+ ///
+ /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition.
+ ///
+ /// The inverse of this matrix.
+ public override Matrix Inverse()
+ {
+ var result = (DenseMatrix)Factors.Clone();
+ Control.LinearAlgebraProvider.LUInverseFactored(result.Data, result.RowCount, Pivots);
+ return result;
+ }
+
+ #region Simple arithmetic of type T
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+
+ ///
+ /// Get value of type T equal to minus one
+ ///
+ /// One value
+ protected sealed override Complex MinusOneValueT
+ {
+ get { return -Complex.One; }
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs
new file mode 100644
index 00000000..dec02790
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs
@@ -0,0 +1,203 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+
+ ///
+ /// A class which encapsulates the functionality of the QR decomposition.
+ /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal matrix
+ /// (its columns are orthogonal unit vectors meaning QTQ = I) and R is an upper triangular matrix
+ /// (also called right triangular matrix).
+ ///
+ ///
+ /// The computation of the QR decomposition is done at construction time by Householder transformation.
+ ///
+ public class DenseQR : QR
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// QR factorization when the constructor is called and cache it's factorization.
+ ///
+ /// The matrix to factor.
+ /// If is null.
+ /// If row count is less then column count
+ public DenseQR(DenseMatrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount < matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ MatrixR = matrix.Clone();
+ MatrixQ = new DenseMatrix(matrix.RowCount);
+ Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data);
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A QR factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // The solution X should have the same number of columns as B
+ if (input.ColumnCount != result.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
+ if (MatrixR.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ // The solution X row dimension is equal to the column dimension of A
+ if (MatrixR.ColumnCount != result.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ var dinput = input as DenseMatrix;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment.");
+ }
+
+ var dresult = result as DenseMatrix;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment.");
+ }
+
+ Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data);
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A QR factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Ax=b where A is an m x n matrix
+ // Check that b is a column vector with m entries
+ if (MatrixR.RowCount != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Check that x is a column vector with n entries
+ if (MatrixR.ColumnCount != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var dinput = input as DenseVector;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment.");
+ }
+
+ var dresult = result as DenseVector;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment.");
+ }
+
+ Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data);
+ }
+
+ #region Simple arithmetic of type T
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Returns the absolute value of a specified number.
+ ///
+ /// A number whose absolute is to be found
+ /// Absolute value
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs
new file mode 100644
index 00000000..88fa07c6
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs
@@ -0,0 +1,217 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+
+ ///
+ /// A class which encapsulates the functionality of the singular value decomposition (SVD) for .
+ /// Suppose M is an m-by-n matrix whose entries are real numbers.
+ /// Then there exists a factorization of the form M = UΣVT where:
+ /// - U is an m-by-m unitary matrix;
+ /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
+ /// - VT denotes transpose of V, an n-by-n unitary matrix;
+ /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
+ /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
+ /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.
+ ///
+ ///
+ /// The computation of the singular value decomposition is done at construction time.
+ ///
+ public class DenseSvd : Svd
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// the singular value decomposition when the constructor is called and cache it's decomposition.
+ ///
+ /// The matrix to factor.
+ /// Compute the singular U and VT vectors or not.
+ /// If is null.
+ /// If SVD algorithm failed to converge with matrix .
+ public DenseSvd(DenseMatrix matrix, bool computeVectors)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ ComputeVectors = computeVectors;
+ var nm = Math.Min(matrix.RowCount, matrix.ColumnCount);
+ VectorS = new DenseVector(nm);
+ MatrixU = new DenseMatrix(matrix.RowCount);
+ MatrixVT = new DenseMatrix(matrix.ColumnCount);
+ Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix)matrix.Clone()).Data, matrix.RowCount, matrix.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data);
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A SVD factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (!ComputeVectors)
+ {
+ throw new InvalidOperationException(Resources.SingularVectorsNotComputed);
+ }
+
+ // The solution X should have the same number of columns as B
+ if (input.ColumnCount != result.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
+ if (MatrixU.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ // The solution X row dimension is equal to the column dimension of A
+ if (MatrixVT.ColumnCount != result.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ var dinput = input as DenseMatrix;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment.");
+ }
+
+ var dresult = result as DenseMatrix;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment.");
+ }
+
+ Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, input.ColumnCount, dresult.Data);
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A SVD factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (!ComputeVectors)
+ {
+ throw new InvalidOperationException(Resources.SingularVectorsNotComputed);
+ }
+
+ // Ax=b where A is an m x n matrix
+ // Check that b is a column vector with m entries
+ if (MatrixU.RowCount != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Check that x is a column vector with n entries
+ if (MatrixVT.ColumnCount != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var dinput = input as DenseVector;
+ if (dinput == null)
+ {
+ throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment.");
+ }
+
+ var dresult = result as DenseVector;
+ if (dresult == null)
+ {
+ throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment.");
+ }
+
+ Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data);
+ }
+
+ #region Simple arithmetic of type T
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Returns the absolute value of a specified number.
+ ///
+ /// A number whose absolute is to be found
+ /// Absolute value
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/GramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/GramSchmidt.cs
new file mode 100644
index 00000000..5e0c85ad
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/GramSchmidt.cs
@@ -0,0 +1,325 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+
+ ///
+ /// A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.
+ /// Any complex square matrix A may be decomposed as A = QR where Q is an unitary mxn matrix and R is an nxn upper triangular matrix.
+ ///
+ ///
+ /// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
+ ///
+ public class GramSchmidt : QR
+ {
+ ///
+ /// Initializes a new instance of the class. This object creates an unitary matrix
+ /// using the modified Gram-Schmidt method.
+ ///
+ /// The matrix to factor.
+ /// If is null.
+ /// If row count is less then column count
+ /// If is rank deficient
+ public GramSchmidt(Matrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount < matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ MatrixQ = matrix.Clone();
+ MatrixR = matrix.CreateMatrix(matrix.ColumnCount, matrix.ColumnCount);
+
+ for (var k = 0; k < MatrixQ.ColumnCount; k++)
+ {
+ var norm = MatrixQ.Column(k).Norm(2);
+ if (norm == 0.0)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient);
+ }
+
+ MatrixR.At(k, k, norm);
+ for (var i = 0; i < MatrixQ.RowCount; i++)
+ {
+ MatrixQ.At(i, k, MatrixQ.At(i, k) / norm);
+ }
+
+ for (var j = k + 1; j < MatrixQ.ColumnCount; j++)
+ {
+ var dot = Complex.Zero;
+ for (int i = 0; i < MatrixQ.RowCount; i++)
+ {
+ dot += MatrixQ.Column(k)[i].Conjugate() * MatrixQ.Column(j)[i];
+ }
+
+ MatrixR.At(k, j, dot);
+ for (var i = 0; i < MatrixQ.RowCount; i++)
+ {
+ var value = MatrixQ.At(i, j) - (MatrixQ.At(i, k) * dot);
+ MatrixQ.At(i, j, value);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Gets a value indicating whether the matrix is full rank or not.
+ ///
+ /// true if the matrix is full rank; otherwise false.
+ public override bool IsFullRank
+ {
+ get
+ {
+ return true;
+ }
+ }
+
+ ///
+ /// Gets the determinant of the matrix for which the QR matrix was computed.
+ ///
+ public override double Determinant
+ {
+ get
+ {
+ if (MatrixQ.RowCount != MatrixQ.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare);
+ }
+
+ var det = Complex.One;
+ for (var i = 0; i < MatrixR.ColumnCount; i++)
+ {
+ det *= MatrixR.At(i, i);
+ if (MatrixR.At(i, i).Magnitude.AlmostEqualInDecimalPlaces(0.0, 15))
+ {
+ return 0;
+ }
+ }
+
+ return det.Magnitude;
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A QR factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // The solution X should have the same number of columns as B
+ if (input.ColumnCount != result.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
+ if (MatrixQ.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ // The solution X row dimension is equal to the column dimension of A
+ if (MatrixQ.ColumnCount != result.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ var inputCopy = input.Clone();
+
+ // Compute Y = transpose(Q)*B
+ var column = new Complex[MatrixQ.RowCount];
+ for (var j = 0; j < input.ColumnCount; j++)
+ {
+ for (var k = 0; k < MatrixQ.RowCount; k++)
+ {
+ column[k] = inputCopy.At(k, j);
+ }
+
+ for (var i = 0; i < MatrixQ.ColumnCount; i++)
+ {
+ var s = Complex.Zero;
+ for (var k = 0; k < MatrixQ.RowCount; k++)
+ {
+ s += MatrixQ.At(k, i).Conjugate() * column[k];
+ }
+
+ inputCopy.At(i, j, s);
+ }
+ }
+
+ // Solve R*X = Y;
+ for (var k = MatrixQ.ColumnCount - 1; k >= 0; k--)
+ {
+ for (var j = 0; j < input.ColumnCount; j++)
+ {
+ inputCopy.At(k, j, inputCopy.At(k, j) / MatrixR.At(k, k));
+ }
+
+ for (var i = 0; i < k; i++)
+ {
+ for (var j = 0; j < input.ColumnCount; j++)
+ {
+ inputCopy.At(i, j, inputCopy.At(i, j) - (inputCopy.At(k, j) * MatrixR.At(i, k)));
+ }
+ }
+ }
+
+ for (var i = 0; i < MatrixR.ColumnCount; i++)
+ {
+ for (var j = 0; j < input.ColumnCount; j++)
+ {
+ result.At(i, j, inputCopy.At(i, j));
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A QR factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Ax=b where A is an m x n matrix
+ // Check that b is a column vector with m entries
+ if (MatrixQ.RowCount != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Check that x is a column vector with n entries
+ if (MatrixQ.ColumnCount != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var inputCopy = input.Clone();
+
+ // Compute Y = transpose(Q)*B
+ var column = new Complex[MatrixQ.RowCount];
+ for (var k = 0; k < MatrixQ.RowCount; k++)
+ {
+ column[k] = inputCopy[k];
+ }
+
+ for (var i = 0; i < MatrixQ.ColumnCount; i++)
+ {
+ var s = Complex.Zero;
+ for (var k = 0; k < MatrixQ.RowCount; k++)
+ {
+ s += MatrixQ.At(k, i).Conjugate() * column[k];
+ }
+
+ inputCopy[i] = s;
+ }
+
+ // Solve R*X = Y;
+ for (var k = MatrixQ.ColumnCount - 1; k >= 0; k--)
+ {
+ inputCopy[k] /= MatrixR.At(k, k);
+ for (var i = 0; i < k; i++)
+ {
+ inputCopy[i] -= inputCopy[k] * MatrixR.At(i, k);
+ }
+ }
+
+ for (var i = 0; i < MatrixR.ColumnCount; i++)
+ {
+ result[i] = inputCopy[i];
+ }
+ }
+
+ #region Simple arithmetic of type T
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Returns the absolute value of a specified number.
+ ///
+ /// A number whose absolute is to be found
+ /// Absolute value
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs
new file mode 100644
index 00000000..1417c631
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs
@@ -0,0 +1,268 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+
+ ///
+ /// A class which encapsulates the functionality of a Cholesky factorization for user matrices.
+ /// For a symmetric, positive definite matrix A, the Cholesky factorization
+ /// is an lower triangular matrix L so that A = L*L'.
+ ///
+ ///
+ /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
+ /// or positive definite, the constructor will throw an exception.
+ ///
+ public class UserCholesky : Cholesky
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// Cholesky factorization when the constructor is called and cache it's factorization.
+ ///
+ /// The matrix to factor.
+ /// If is null.
+ /// If is not a square matrix.
+ /// If is not positive definite.
+ public UserCholesky(Matrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare);
+ }
+
+ // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
+ CholeskyFactor = matrix.Clone();
+ for (var i = 0; i < CholeskyFactor.RowCount; i++)
+ {
+ var d = Complex.Zero;
+ for (var j = 0; j < i; j++)
+ {
+ var s = Complex.Zero;
+ for (var k = 0; k < j; k++)
+ {
+ s += CholeskyFactor.At(i, k) * CholeskyFactor.At(j, k).Conjugate();
+ }
+
+ s = (matrix.At(i, j) - s) / CholeskyFactor.At(j, j);
+ CholeskyFactor.At(i, j, s);
+ d += s * s.Conjugate();
+ }
+
+ d = matrix.At(i, i) - d;
+ if (d.Real <= 0)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
+ }
+
+ CholeskyFactor.At(i, i, d.SquareRoot());
+ for (var k = i + 1; k < CholeskyFactor.RowCount; k++)
+ {
+ CholeskyFactor.At(i, k, 0.0);
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A Cholesky factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (result.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ if (result.ColumnCount != input.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ if (input.RowCount != CholeskyFactor.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ input.CopyTo(result);
+ var order = CholeskyFactor.RowCount;
+
+ for (var c = 0; c < result.ColumnCount; c++)
+ {
+ // Solve L*Y = B;
+ Complex sum;
+ for (var i = 0; i < order; i++)
+ {
+ sum = result.At(i, c);
+ for (var k = i - 1; k >= 0; k--)
+ {
+ sum -= CholeskyFactor.At(i, k) * result.At(k, c);
+ }
+
+ result.At(i, c, sum / CholeskyFactor.At(i, i));
+ }
+
+ // Solve L'*X = Y;
+ for (var i = order - 1; i >= 0; i--)
+ {
+ sum = result.At(i, c);
+ for (var k = i + 1; k < order; k++)
+ {
+ sum -= CholeskyFactor.At(k, i).Conjugate() * result.At(k, c);
+ }
+
+ result.At(i, c, sum / CholeskyFactor.At(i, i));
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A Cholesky factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (input.Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ if (input.Count != CholeskyFactor.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ input.CopyTo(result);
+ var order = CholeskyFactor.RowCount;
+
+ // Solve L*Y = B;
+ Complex sum;
+ for (var i = 0; i < order; i++)
+ {
+ sum = result[i];
+ for (var k = i - 1; k >= 0; k--)
+ {
+ sum -= CholeskyFactor.At(i, k) * result[k];
+ }
+
+ result[i] = sum / CholeskyFactor.At(i, i);
+ }
+
+ // Solve L'*X = Y;
+ for (var i = order - 1; i >= 0; i--)
+ {
+ sum = result[i];
+ for (var k = i + 1; k < order; k++)
+ {
+ sum -= CholeskyFactor.At(k, i).Conjugate() * result[k];
+ }
+
+ result[i] = sum / CholeskyFactor.At(i, i);
+ }
+ }
+
+ #region Simple arithmetic of type T
+ ///
+ /// Add two values T+T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of addition
+ protected sealed override Complex AddT(Complex val1, Complex val2)
+ {
+ return val1 + val2;
+ }
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Returns the natural (base e) logarithm of a specified number.
+ ///
+ /// A number whose logarithm is to be found
+ /// Natural (base e) logarithm
+ protected sealed override Complex LogT(Complex val1)
+ {
+ return val1.NaturalLogarithm();
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs
new file mode 100644
index 00000000..a0768b70
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs
@@ -0,0 +1,335 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+
+ ///
+ /// A class which encapsulates the functionality of an LU factorization.
+ /// For a matrix A, the LU factorization is a pair of lower triangular matrix L and
+ /// upper triangular matrix U so that A = L*U.
+ ///
+ ///
+ /// The computation of the LU factorization is done at construction time.
+ ///
+ public class UserLU : LU
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// LU factorization when the constructor is called and cache it's factorization.
+ ///
+ /// The matrix to factor.
+ /// If is null.
+ /// If is not a square matrix.
+ public UserLU(Matrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare);
+ }
+
+ // Create an array for the pivot indices.
+ var order = matrix.RowCount;
+ Factors = matrix.Clone();
+ Pivots = new int[order];
+
+ // Initialize the pivot matrix to the identity permutation.
+ for (var i = 0; i < order; i++)
+ {
+ Pivots[i] = i;
+ }
+
+ var vectorLUcolj = new Complex[order];
+ for (var j = 0; j < order; j++)
+ {
+ // Make a copy of the j-th column to localize references.
+ for (var i = 0; i < order; i++)
+ {
+ vectorLUcolj[i] = Factors.At(i, j);
+ }
+
+ // Apply previous transformations.
+ for (var i = 0; i < order; i++)
+ {
+ var kmax = Math.Min(i, j);
+ var s = Complex.Zero;
+ for (var k = 0; k < kmax; k++)
+ {
+ s += Factors.At(i, k) * vectorLUcolj[k];
+ }
+
+ vectorLUcolj[i] -= s;
+ Factors.At(i, j, vectorLUcolj[i]);
+ }
+
+ // Find pivot and exchange if necessary.
+ var p = j;
+ for (var i = j + 1; i < order; i++)
+ {
+ if (vectorLUcolj[i].Magnitude > vectorLUcolj[p].Magnitude)
+ {
+ p = i;
+ }
+ }
+
+ if (p != j)
+ {
+ for (var k = 0; k < order; k++)
+ {
+ var temp = Factors.At(p, k);
+ Factors.At(p, k, Factors.At(j, k));
+ Factors.At(j, k, temp);
+ }
+
+ Pivots[j] = p;
+ }
+
+ // Compute multipliers.
+ if (j < order & Factors.At(j, j) != 0.0)
+ {
+ for (var i = j + 1; i < order; i++)
+ {
+ Factors.At(i, j, (Factors.At(i, j) / Factors.At(j, j)));
+ }
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A LU factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (result.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ if (result.ColumnCount != input.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ if (input.RowCount != Factors.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ // Copy the contents of input to result.
+ input.CopyTo(result);
+ for (var i = 0; i < Pivots.Length; i++)
+ {
+ if (Pivots[i] == i)
+ {
+ continue;
+ }
+
+ var p = Pivots[i];
+ for (var j = 0; j < result.ColumnCount; j++)
+ {
+ var temp = result.At(p, j);
+ result.At(p, j, result.At(i, j));
+ result.At(i, j, temp);
+ }
+ }
+
+ var order = Factors.RowCount;
+
+ // Solve L*Y = P*B
+ for (var k = 0; k < order; k++)
+ {
+ for (var i = k + 1; i < order; i++)
+ {
+ for (var j = 0; j < result.ColumnCount; j++)
+ {
+ var temp = result.At(k, j) * Factors.At(i, k);
+ result.At(i, j, result.At(i, j) - temp);
+ }
+ }
+ }
+
+ // Solve U*X = Y;
+ for (var k = order - 1; k >= 0; k--)
+ {
+ for (var j = 0; j < result.ColumnCount; j++)
+ {
+ result.At(k, j, (result.At(k, j) / Factors.At(k, k)));
+ }
+
+ for (var i = 0; i < k; i++)
+ {
+ for (var j = 0; j < result.ColumnCount; j++)
+ {
+ var temp = result.At(k, j) * Factors.At(i, k);
+ result.At(i, j, result.At(i, j) - temp);
+ }
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A LU factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Check for proper dimensions.
+ if (input.Count != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ if (input.Count != Factors.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ // Copy the contents of input to result.
+ input.CopyTo(result);
+ for (var i = 0; i < Pivots.Length; i++)
+ {
+ if (Pivots[i] == i)
+ {
+ continue;
+ }
+
+ var p = Pivots[i];
+ var temp = result[p];
+ result[p] = result[i];
+ result[i] = temp;
+ }
+
+ var order = Factors.RowCount;
+
+ // Solve L*Y = P*B
+ for (var k = 0; k < order; k++)
+ {
+ for (var i = k + 1; i < order; i++)
+ {
+ result[i] -= result[k] * Factors.At(i, k);
+ }
+ }
+
+ // Solve U*X = Y;
+ for (var k = order - 1; k >= 0; k--)
+ {
+ result[k] /= Factors.At(k, k);
+ for (var i = 0; i < k; i++)
+ {
+ result[i] -= result[k] * Factors.At(i, k);
+ }
+ }
+ }
+
+ ///
+ /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition.
+ ///
+ /// The inverse of this matrix.
+ public override Matrix Inverse()
+ {
+ var order = Factors.RowCount;
+ var inverse = Factors.CreateMatrix(order, order);
+ for (var i = 0; i < order; i++)
+ {
+ inverse.At(i, i, 1.0);
+ }
+
+ return Solve(inverse);
+ }
+
+ #region Simple arithmetic of type T
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+
+ ///
+ /// Get value of type T equal to minus one
+ ///
+ /// One value
+ protected sealed override Complex MinusOneValueT
+ {
+ get { return -Complex.One; }
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs
new file mode 100644
index 00000000..293fb4f5
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs
@@ -0,0 +1,366 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Linq;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+
+ ///
+ /// A class which encapsulates the functionality of the QR decomposition.
+ /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal matrix
+ /// (its columns are orthogonal unit vectors meaning QTQ = I) and R is an upper triangular matrix
+ /// (also called right triangular matrix).
+ ///
+ ///
+ /// The computation of the QR decomposition is done at construction time by Householder transformation.
+ ///
+ public class UserQR : QR
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// QR factorization when the constructor is called and cache it's factorization.
+ ///
+ /// The matrix to factor.
+ /// If is null.
+ public UserQR(Matrix matrix)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount < matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ MatrixR = matrix.Clone();
+ MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount);
+
+ for (var i = 0; i < matrix.RowCount; i++)
+ {
+ MatrixQ.At(i, i, 1.0);
+ }
+
+ var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
+ var u = new Complex[minmn][];
+ for (var i = 0; i < minmn; i++)
+ {
+ u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
+ ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
+ }
+
+ for (var i = minmn - 1; i >= 0; i--)
+ {
+ ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
+ }
+ }
+
+ ///
+ /// Generate column from initial matrix to work array
+ ///
+ /// Initial matrix
+ /// The firts row
+ /// The last row
+ /// Column index
+ /// Generated vector
+ private static Complex[] GenerateColumn(Matrix a, int rowStart, int rowEnd, int column)
+ {
+ var ru = rowEnd - rowStart + 1;
+ var u = new Complex[ru];
+
+ for (var i = rowStart; i <= rowEnd; i++)
+ {
+ u[i - rowStart] = a.At(i, column);
+ a.At(i, column, 0.0);
+ }
+
+ var norm = u.Aggregate(Complex.Zero, (current, t) => current + (t.Magnitude * t.Magnitude));
+ norm = norm.SquareRoot();
+
+ if (rowStart == rowEnd || norm.Magnitude == 0)
+ {
+ a.At(rowStart, column, -u[0]);
+ u[0] = Math.Sqrt(2.0);
+ return u;
+ }
+
+ if (u[0].Magnitude != 0.0)
+ {
+ norm = norm.Magnitude * (u[0] / u[0].Magnitude);
+ }
+
+ a.At(rowStart, column, -norm);
+
+ for (var i = 0; i < ru; i++)
+ {
+ u[i] /= norm;
+ }
+
+ u[0] += 1.0;
+
+ var s = (1.0 / u[0]).SquareRoot();
+ for (var i = 0; i < ru; i++)
+ {
+ u[i] = u[i].Conjugate() * s;
+ }
+
+ return u;
+ }
+
+ ///
+ /// Compute matrix Q or R
+ ///
+ /// H vectors values
+ /// Matrix to calculate
+ /// Starting row index
+ /// Ending row index
+ /// Starting column index
+ /// Ending column index
+ private static void ComputeQR(Complex[] u, Matrix a, int rowStart, int rowEnd, int columnStart, int columnEnd)
+ {
+ if (rowEnd < rowStart || columnEnd < columnStart)
+ {
+ return;
+ }
+
+ var v = new Complex[columnEnd - columnStart + 1];
+ for (var j = columnStart; j <= columnEnd; j++)
+ {
+ v[j - columnStart] = 0.0;
+ }
+
+ for (var i = rowStart; i <= rowEnd; i++)
+ {
+ for (var j = columnStart; j <= columnEnd; j++)
+ {
+ v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
+ }
+ }
+
+ for (var i = rowStart; i <= rowEnd; i++)
+ {
+ for (var j = columnStart; j <= columnEnd; j++)
+ {
+ a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * v[j - columnStart]));
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A QR factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // The solution X should have the same number of columns as B
+ if (input.ColumnCount != result.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
+ if (MatrixR.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ // The solution X row dimension is equal to the column dimension of A
+ if (MatrixR.ColumnCount != result.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ var inputCopy = input.Clone();
+
+ // Compute Y = transpose(Q)*B
+ var column = new Complex[MatrixR.RowCount];
+ for (var j = 0; j < input.ColumnCount; j++)
+ {
+ for (var k = 0; k < MatrixR.RowCount; k++)
+ {
+ column[k] = inputCopy.At(k, j);
+ }
+
+ for (var i = 0; i < MatrixR.RowCount; i++)
+ {
+ var s = Complex.Zero;
+ for (var k = 0; k < MatrixR.RowCount; k++)
+ {
+ s += MatrixQ.At(k, i).Conjugate() * column[k];
+ }
+
+ inputCopy.At(i, j, s);
+ }
+ }
+
+ // Solve R*X = Y;
+ for (var k = MatrixR.ColumnCount - 1; k >= 0; k--)
+ {
+ for (var j = 0; j < input.ColumnCount; j++)
+ {
+ inputCopy.At(k, j, inputCopy.At(k, j) / MatrixR.At(k, k));
+ }
+
+ for (var i = 0; i < k; i++)
+ {
+ for (var j = 0; j < input.ColumnCount; j++)
+ {
+ inputCopy.At(i, j, inputCopy.At(i, j) - (inputCopy.At(k, j) * MatrixR.At(i, k)));
+ }
+ }
+ }
+
+ for (var i = 0; i < MatrixR.ColumnCount; i++)
+ {
+ for (var j = 0; j < inputCopy.ColumnCount; j++)
+ {
+ result.At(i, j, inputCopy.At(i, j));
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A QR factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ // Ax=b where A is an m x n matrix
+ // Check that b is a column vector with m entries
+ if (MatrixR.RowCount != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Check that x is a column vector with n entries
+ if (MatrixR.ColumnCount != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var inputCopy = input.Clone();
+
+ // Compute Y = transpose(Q)*B
+ var column = new Complex[MatrixR.RowCount];
+ for (var k = 0; k < MatrixR.RowCount; k++)
+ {
+ column[k] = inputCopy[k];
+ }
+
+ for (var i = 0; i < MatrixR.RowCount; i++)
+ {
+ var s = Complex.Zero;
+ for (var k = 0; k < MatrixR.RowCount; k++)
+ {
+ s += MatrixQ.At(k, i).Conjugate() * column[k];
+ }
+
+ inputCopy[i] = s;
+ }
+
+ // Solve R*X = Y;
+ for (var k = MatrixR.ColumnCount - 1; k >= 0; k--)
+ {
+ inputCopy[k] /= MatrixR.At(k, k);
+ for (var i = 0; i < k; i++)
+ {
+ inputCopy[i] -= inputCopy[k] * MatrixR.At(i, k);
+ }
+ }
+
+ for (var i = 0; i < MatrixR.ColumnCount; i++)
+ {
+ result[i] = inputCopy[i];
+ }
+ }
+
+ #region Simple arithmetic of type T
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Returns the absolute value of a specified number.
+ ///
+ /// A number whose absolute is to be found
+ /// Absolute value
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs
new file mode 100644
index 00000000..0b2ab5fa
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserSvd.cs
@@ -0,0 +1,974 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Factorization;
+ using Properties;
+
+ ///
+ /// A class which encapsulates the functionality of the singular value decomposition (SVD) for .
+ /// Suppose M is an m-by-n matrix whose entries are real numbers.
+ /// Then there exists a factorization of the form M = UΣVT where:
+ /// - U is an m-by-m unitary matrix;
+ /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
+ /// - VT denotes transpose of V, an n-by-n unitary matrix;
+ /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
+ /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
+ /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.
+ ///
+ ///
+ /// The computation of the singular value decomposition is done at construction time.
+ ///
+ public class UserSvd : Svd
+ {
+ ///
+ /// Initializes a new instance of the class. This object will compute the
+ /// the singular value decomposition when the constructor is called and cache it's decomposition.
+ ///
+ /// The matrix to factor.
+ /// Compute the singular U and VT vectors or not.
+ /// If is null.
+ /// If SVD algorithm failed to converge with matrix .
+ public UserSvd(Matrix matrix, bool computeVectors)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ ComputeVectors = computeVectors;
+ var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount);
+ var matrixCopy = matrix.Clone();
+
+ VectorS = matrixCopy.CreateVector(nm);
+ MatrixU = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount);
+ MatrixVT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount);
+
+ const int Maxiter = 1000;
+ var e = new Complex[matrixCopy.ColumnCount];
+ var work = new Complex[matrixCopy.RowCount];
+
+ int i, j;
+ int l, lp1;
+ var cs = 0.0;
+ var sn = 0.0;
+ Complex t;
+
+ var ncu = matrixCopy.RowCount;
+
+ // Reduce matrixCopy to bidiagonal form, storing the diagonal elements
+ // In s and the super-diagonal elements in e.
+ var nct = Math.Min(matrixCopy.RowCount - 1, matrixCopy.ColumnCount);
+ var nrt = Math.Max(0, Math.Min(matrixCopy.ColumnCount - 2, matrixCopy.RowCount));
+ var lu = Math.Max(nct, nrt);
+ for (l = 0; l < lu; l++)
+ {
+ lp1 = l + 1;
+ if (l < nct)
+ {
+ // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l].
+ VectorS[l] = Cnrm2Column(matrixCopy, matrixCopy.RowCount, l, l);
+ if (VectorS[l].Magnitude != 0.0)
+ {
+ if (matrixCopy.At(l, l).Magnitude != 0.0)
+ {
+ VectorS[l] = Csign(VectorS[l], matrixCopy.At(l, l));
+ }
+
+ CscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0 / VectorS[l]);
+ matrixCopy.At(l, l, (Complex.One + matrixCopy.At(l, l)));
+ }
+
+ VectorS[l] = -VectorS[l];
+ }
+
+ for (j = lp1; j < matrixCopy.ColumnCount; j++)
+ {
+ if (l < nct)
+ {
+ if (VectorS[l].Magnitude != 0.0)
+ {
+ // Apply the transformation.
+ t = -Cdotc(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l);
+ if (t != Complex.Zero)
+ {
+ for (var ii = l; ii < matrixCopy.RowCount; ii++)
+ {
+ matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l)));
+ }
+ }
+ }
+ }
+
+ // Place the l-th row of matrixCopy into e for the
+ // Subsequent calculation of the row transformation.
+ e[j] = matrixCopy.At(l, j).Conjugate();
+ }
+
+ if (ComputeVectors && l < nct)
+ {
+ // Place the transformation in u for subsequent back multiplication.
+ for (i = l; i < matrixCopy.RowCount; i++)
+ {
+ MatrixU.At(i, l, matrixCopy.At(i, l));
+ }
+ }
+
+ if (l >= nrt)
+ {
+ continue;
+ }
+
+ // Compute the l-th row transformation and place the l-th super-diagonal in e(l).
+ var enorm = Cnrm2Vector(e, lp1);
+ e[l] = enorm;
+ if (e[l].Magnitude != 0.0)
+ {
+ if (e[lp1].Magnitude != 0.0)
+ {
+ e[l] = Csign(e[l], e[lp1]);
+ }
+
+ CscalVector(e, lp1, 1.0 / e[l]);
+ e[lp1] = Complex.One + e[lp1];
+ }
+
+ e[l] = -e[l].Conjugate();
+ if (lp1 < matrixCopy.RowCount && e[l].Magnitude != 0.0)
+ {
+ // Apply the transformation.
+ for (i = lp1; i < matrixCopy.RowCount; i++)
+ {
+ work[i] = Complex.Zero;
+ }
+
+ for (j = lp1; j < matrixCopy.ColumnCount; j++)
+ {
+ if (e[j] != Complex.Zero)
+ {
+ for (var ii = lp1; ii < matrixCopy.RowCount; ii++)
+ {
+ work[ii] += e[j] * matrixCopy.At(ii, j);
+ }
+ }
+ }
+
+ for (j = lp1; j < matrixCopy.ColumnCount; j++)
+ {
+ var ww = (-e[j] / e[lp1]).Conjugate();
+ if (ww != Complex.Zero)
+ {
+ for (var ii = lp1; ii < matrixCopy.RowCount; ii++)
+ {
+ matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii]));
+ }
+ }
+ }
+ }
+
+ if (ComputeVectors)
+ {
+ // Place the transformation in v for subsequent back multiplication.
+ for (i = lp1; i < matrixCopy.ColumnCount; i++)
+ {
+ MatrixVT.At(i, l, e[i]);
+ }
+ }
+ }
+
+ // Set up the final bidiagonal matrixCopy or order m.
+ var m = Math.Min(matrixCopy.ColumnCount, matrixCopy.RowCount + 1);
+ var nctp1 = nct + 1;
+ var nrtp1 = nrt + 1;
+ if (nct < matrixCopy.ColumnCount)
+ {
+ VectorS[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1));
+ }
+
+ if (matrixCopy.RowCount < m)
+ {
+ VectorS[m - 1] = Complex.Zero;
+ }
+
+ if (nrtp1 < m)
+ {
+ e[nrtp1 - 1] = matrixCopy.At((nrtp1 - 1), (m - 1));
+ }
+
+ e[m - 1] = Complex.Zero;
+
+ // If required, generate u.
+ if (ComputeVectors)
+ {
+ for (j = nctp1 - 1; j < ncu; j++)
+ {
+ for (i = 0; i < matrixCopy.RowCount; i++)
+ {
+ MatrixU.At(i, j, Complex.Zero);
+ }
+
+ MatrixU.At(j, j, Complex.One);
+ }
+
+ for (l = nct - 1; l >= 0; l--)
+ {
+ if (VectorS[l].Magnitude != 0.0)
+ {
+ for (j = l + 1; j < ncu; j++)
+ {
+ t = -Cdotc(MatrixU, matrixCopy.RowCount, l, j, l) / MatrixU.At(l, l);
+ if (t != Complex.Zero)
+ {
+ for (var ii = l; ii < matrixCopy.RowCount; ii++)
+ {
+ MatrixU.At(ii, j, MatrixU.At(ii, j) + (t * MatrixU.At(ii, l)));
+ }
+ }
+ }
+
+ CscalColumn(MatrixU, matrixCopy.RowCount, l, l, -1.0);
+ MatrixU.At(l, l, Complex.One + MatrixU.At(l, l));
+ for (i = 0; i < l; i++)
+ {
+ MatrixU.At(i, l, Complex.Zero);
+ }
+ }
+ else
+ {
+ for (i = 0; i < matrixCopy.RowCount; i++)
+ {
+ MatrixU.At(i, l, Complex.Zero);
+ }
+
+ MatrixU.At(l, l, Complex.One);
+ }
+ }
+ }
+
+ // If it is required, generate v.
+ if (ComputeVectors)
+ {
+ for (l = matrixCopy.ColumnCount - 1; l >= 0; l--)
+ {
+ lp1 = l + 1;
+ if (l < nrt)
+ {
+ if (e[l].Magnitude != 0.0)
+ {
+ for (j = lp1; j < matrixCopy.ColumnCount; j++)
+ {
+ t = -Cdotc(MatrixVT, matrixCopy.ColumnCount, l, j, lp1) / MatrixVT.At(lp1, l);
+ if (t != Complex.Zero)
+ {
+ for (var ii = l; ii < matrixCopy.ColumnCount; ii++)
+ {
+ MatrixVT.At(ii, j, MatrixVT.At(ii, j) + (t * MatrixVT.At(ii, l)));
+ }
+ }
+ }
+ }
+ }
+
+ for (i = 0; i < matrixCopy.ColumnCount; i++)
+ {
+ MatrixVT.At(i, l, Complex.Zero);
+ }
+
+ MatrixVT.At(l, l, Complex.One);
+ }
+ }
+
+ // Transform s and e so that they are real .
+ for (i = 0; i < m; i++)
+ {
+ Complex r;
+ if (VectorS[i].Magnitude != 0.0)
+ {
+ t = VectorS[i].Magnitude;
+ r = VectorS[i] / t;
+ VectorS[i] = t;
+ if (i < m - 1)
+ {
+ e[i] = e[i] / r;
+ }
+
+ if (ComputeVectors)
+ {
+ CscalColumn(MatrixU, matrixCopy.RowCount, i, 0, r);
+ }
+ }
+
+ // Exit
+ if (i == m - 1)
+ {
+ break;
+ }
+
+ if (e[i].Magnitude != 0.0)
+ {
+ t = e[i].Magnitude;
+ r = t / e[i];
+ e[i] = t;
+ VectorS[i + 1] = VectorS[i + 1] * r;
+ if (ComputeVectors)
+ {
+ CscalColumn(MatrixVT, matrixCopy.ColumnCount, i + 1, 0, r);
+ }
+ }
+ }
+
+ // Main iteration loop for the singular values.
+ var mn = m;
+ var iter = 0;
+
+ while (m > 0)
+ {
+ // Quit if all the singular values have been found. If too many iterations have been performed,
+ // throw exception that Convergence Failed
+ if (iter >= Maxiter)
+ {
+ throw new ArgumentException(Resources.ConvergenceFailed);
+ }
+
+ // This section of the program inspects for negligible elements in the s and e arrays. On
+ // completion the variables kase and l are set as follows.
+ // Kase = 1 if VectorS[m] and e[l-1] are negligible and l < m
+ // Kase = 2 if VectorS[l] is negligible and l < m
+ // Kase = 3 if e[l-1] is negligible, l < m, and VectorS[l, ..., VectorS[m] are not negligible (qr step).
+ // Лase = 4 if e[m-1] is negligible (convergence).
+ double ztest;
+ double test;
+ for (l = m - 2; l >= 0; l--)
+ {
+ test = VectorS[l].Magnitude + VectorS[l + 1].Magnitude;
+ ztest = test + e[l].Magnitude;
+ if (ztest.AlmostEqualInDecimalPlaces(test, 15))
+ {
+ e[l] = Complex.Zero;
+ break;
+ }
+ }
+
+ int kase;
+ if (l == m - 2)
+ {
+ kase = 4;
+ }
+ else
+ {
+ int ls;
+ for (ls = m - 1; ls > l; ls--)
+ {
+ test = 0.0;
+ if (ls != m - 1)
+ {
+ test = test + e[ls].Magnitude;
+ }
+
+ if (ls != l + 1)
+ {
+ test = test + e[ls - 1].Magnitude;
+ }
+
+ ztest = test + VectorS[ls].Magnitude;
+ if (ztest.AlmostEqualInDecimalPlaces(test, 15))
+ {
+ VectorS[ls] = Complex.Zero;
+ break;
+ }
+ }
+
+ if (ls == l)
+ {
+ kase = 3;
+ }
+ else if (ls == m - 1)
+ {
+ kase = 1;
+ }
+ else
+ {
+ kase = 2;
+ l = ls;
+ }
+ }
+
+ l = l + 1;
+
+ // Perform the task indicated by kase.
+ int k;
+ double f;
+ switch (kase)
+ {
+ // Deflate negligible VectorS[m].
+ case 1:
+ f = e[m - 2].Real;
+ e[m - 2] = Complex.Zero;
+ double t1;
+ for (var kk = l; kk < m - 1; kk++)
+ {
+ k = m - 2 - kk + l;
+ t1 = VectorS[k].Real;
+ Srotg(ref t1, ref f, ref cs, ref sn);
+ VectorS[k] = t1;
+ if (k != l)
+ {
+ f = -sn * e[k - 1].Real;
+ e[k - 1] = cs * e[k - 1];
+ }
+
+ if (ComputeVectors)
+ {
+ Csrot(MatrixVT, matrixCopy.ColumnCount, k, m - 1, cs, sn);
+ }
+ }
+
+ break;
+
+ // Split at negligible VectorS[l].
+ case 2:
+ f = e[l - 1].Real;
+ e[l - 1] = Complex.Zero;
+ for (k = l; k < m; k++)
+ {
+ t1 = VectorS[k].Real;
+ Srotg(ref t1, ref f, ref cs, ref sn);
+ VectorS[k] = t1;
+ f = -sn * e[k].Real;
+ e[k] = cs * e[k];
+ if (ComputeVectors)
+ {
+ Csrot(MatrixU, matrixCopy.RowCount, k, l - 1, cs, sn);
+ }
+ }
+
+ break;
+
+ // Perform one qr step.
+ case 3:
+ // Calculate the shift.
+ var scale = 0.0;
+ scale = Math.Max(scale, VectorS[m - 1].Magnitude);
+ scale = Math.Max(scale, VectorS[m - 2].Magnitude);
+ scale = Math.Max(scale, e[m - 2].Magnitude);
+ scale = Math.Max(scale, VectorS[l].Magnitude);
+ scale = Math.Max(scale, e[l].Magnitude);
+ var sm = VectorS[m - 1].Real / scale;
+ var smm1 = VectorS[m - 2].Real / scale;
+ var emm1 = e[m - 2].Real / scale;
+ var sl = VectorS[l].Real / scale;
+ var el = e[l].Real / scale;
+ var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0;
+ var c = (sm * emm1) * (sm * emm1);
+ var shift = 0.0;
+
+ if (b != 0.0 || c != 0.0)
+ {
+ shift = Math.Sqrt((b * b) + c);
+ if (b < 0.0)
+ {
+ shift = -shift;
+ }
+
+ shift = c / (b + shift);
+ }
+
+ f = ((sl + sm) * (sl - sm)) + shift;
+ var g = sl * el;
+
+ // Chase zeros.
+ for (k = l; k < m - 1; k++)
+ {
+ Srotg(ref f, ref g, ref cs, ref sn);
+ if (k != l)
+ {
+ e[k - 1] = f;
+ }
+
+ f = (cs * VectorS[k].Real) + (sn * e[k].Real);
+ e[k] = (cs * e[k]) - (sn * VectorS[k]);
+ g = sn * VectorS[k + 1].Real;
+ VectorS[k + 1] = cs * VectorS[k + 1];
+ if (ComputeVectors)
+ {
+ Csrot(MatrixVT, matrixCopy.ColumnCount, k, k + 1, cs, sn);
+ }
+
+ Srotg(ref f, ref g, ref cs, ref sn);
+ VectorS[k] = f;
+ f = (cs * e[k].Real) + (sn * VectorS[k + 1].Real);
+ VectorS[k + 1] = (-sn * e[k]) + (cs * VectorS[k + 1]);
+ g = sn * e[k + 1].Real;
+ e[k + 1] = cs * e[k + 1];
+ if (ComputeVectors && k < matrixCopy.RowCount)
+ {
+ Csrot(MatrixU, matrixCopy.RowCount, k, k + 1, cs, sn);
+ }
+ }
+
+ e[m - 2] = f;
+ iter = iter + 1;
+ break;
+
+ // Convergence.
+ case 4:
+ // Make the singular value positive
+ if (VectorS[l].Real < 0.0)
+ {
+ VectorS[l] = -VectorS[l];
+ if (ComputeVectors)
+ {
+ CscalColumn(MatrixVT, matrixCopy.ColumnCount, l, 0, -1.0);
+ }
+ }
+
+ // Order the singular value.
+ while (l != mn - 1)
+ {
+ if (VectorS[l].Real >= VectorS[l + 1].Real)
+ {
+ break;
+ }
+
+ t = VectorS[l];
+ VectorS[l] = VectorS[l + 1];
+ VectorS[l + 1] = t;
+ if (ComputeVectors && l < matrixCopy.ColumnCount)
+ {
+ Swap(MatrixVT, matrixCopy.ColumnCount, l, l + 1);
+ }
+
+ if (ComputeVectors && l < matrixCopy.RowCount)
+ {
+ Swap(MatrixU, matrixCopy.RowCount, l, l + 1);
+ }
+
+ l = l + 1;
+ }
+
+ iter = 0;
+ m = m - 1;
+ break;
+ }
+ }
+
+ if (ComputeVectors)
+ {
+ MatrixVT = MatrixVT.ConjugateTranspose();
+ }
+
+ // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses
+ // a singular vector of length mRows+1 when mRows < mColumns. The last element is not used and needs to be removed.
+ // we should port lapack's svd routine to remove this problem.
+ if (matrixCopy.RowCount < matrixCopy.ColumnCount)
+ {
+ nm--;
+ var tmp = matrixCopy.CreateVector(nm);
+ for (i = 0; i < nm; i++)
+ {
+ tmp[i] = VectorS[i];
+ }
+
+ VectorS = tmp;
+ }
+ }
+
+ ///
+ /// Calculates absolute value of multiplied on signum function of
+ ///
+ /// Complex value z1
+ /// Complex value z2
+ /// Result multiplication of signum function and absolute value
+ private static Complex Csign(Complex z1, Complex z2)
+ {
+ return z1.Magnitude * (z2 / z2.Magnitude);
+ }
+
+ ///
+ /// Interchanges two vectors and
+ ///
+ /// Source matrix
+ /// The number of rows in
+ /// Column A index to swap
+ /// Column B index to swap
+ private static void Swap(Matrix a, int rowCount, int columnA, int columnB)
+ {
+ for (var i = 0; i < rowCount; i++)
+ {
+ var z = a.At(i, columnA);
+ a.At(i, columnA, a.At(i, columnB));
+ a.At(i, columnB, z);
+ }
+ }
+
+ ///
+ /// Scale column by starting from row
+ ///
+ /// Source matrix
+ /// The number of rows in
+ /// Column to scale
+ /// Row to scale from
+ /// Scale value
+ private static void CscalColumn(Matrix a, int rowCount, int column, int rowStart, Complex z)
+ {
+ for (var i = rowStart; i < rowCount; i++)
+ {
+ a.At(i, column, a.At(i, column) * z);
+ }
+ }
+
+ ///
+ /// Scale vector by starting from index
+ ///
+ /// Source vector
+ /// Row to scale from
+ /// Scale value
+ private static void CscalVector(Complex[] a, int start, Complex z)
+ {
+ for (var i = start; i < a.Length; i++)
+ {
+ a[i] = a[i] * z;
+ }
+ }
+
+ ///
+ /// Given the Cartesian coordinates (da, db) of a point p, these fucntion return the parameters da, db, c, and s
+ /// associated with the Givens rotation that zeros the y-coordinate of the point.
+ ///
+ /// Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation
+ /// Provides the y-coordinate of the point p. On exit contains the parameter z associated with the Givens rotation
+ /// Contains the parameter c associated with the Givens rotation
+ /// Contains the parameter s associated with the Givens rotation
+ /// This is equivalent to the DROTG LAPACK routine.
+ private static void Srotg(ref double da, ref double db, ref double c, ref double s)
+ {
+ double r, z;
+
+ var roe = db;
+ var absda = Math.Abs(da);
+ var absdb = Math.Abs(db);
+ if (absda > absdb)
+ {
+ roe = da;
+ }
+
+ var scale = absda + absdb;
+ if (scale == 0.0)
+ {
+ c = 1.0;
+ s = 0.0;
+ r = 0.0;
+ z = 0.0;
+ }
+ else
+ {
+ var sda = da / scale;
+ var sdb = db / scale;
+ r = scale * Math.Sqrt((sda * sda) + (sdb * sdb));
+ if (roe < 0.0)
+ {
+ r = -r;
+ }
+
+ c = da / r;
+ s = db / r;
+ z = 1.0;
+ if (absda > absdb)
+ {
+ z = s;
+ }
+
+ if (absdb >= absda && c != 0.0)
+ {
+ z = 1.0 / c;
+ }
+ }
+
+ da = r;
+ db = z;
+ }
+
+ /// dded
+ /// Calculate Norm 2 of the column in matrix starting from row
+ ///
+ /// Source matrix
+ /// The number of rows in
+ /// Column index
+ /// Start row index
+ /// Norm2 (Euclidean norm) of trhe column
+ private static double Cnrm2Column(Matrix a, int rowCount, int column, int rowStart)
+ {
+ var s = 0.0;
+ for (var i = rowStart; i < rowCount; i++)
+ {
+ s += a.At(i, column).Magnitude * a.At(i, column).Magnitude;
+ }
+
+ return Math.Sqrt(s);
+ }
+
+ ///
+ /// Calculate Norm 2 of the vector starting from index
+ ///
+ /// Source vector
+ /// Start index
+ /// Norm2 (Euclidean norm) of the vector
+ private static double Cnrm2Vector(Complex[] a, int rowStart)
+ {
+ var s = 0.0;
+ for (var i = rowStart; i < a.Length; i++)
+ {
+ s += a[i].Magnitude * a[i].Magnitude;
+ }
+
+ return Math.Sqrt(s);
+ }
+
+ ///
+ /// Calculate dot product of and conjugating the first vector.
+ ///
+ /// Source matrix
+ /// The number of rows in
+ /// Index of column A
+ /// Index of column B
+ /// Starting row index
+ /// Dot product value
+ private static Complex Cdotc(Matrix a, int rowCount, int columnA, int columnB, int rowStart)
+ {
+ var z = Complex.Zero;
+ for (var i = rowStart; i < rowCount; i++)
+ {
+ z += a.At(i, columnA).Conjugate() * a.At(i, columnB);
+ }
+
+ return z;
+ }
+
+ ///
+ /// Performs rotation of points in the plane. Given two vectors x and y ,
+ /// each vector element of these vectors is replaced as follows: x(i) = c*x(i) + s*y(i); y(i) = c*y(i) - s*x(i)
+ ///
+ /// Source matrix
+ /// The number of rows in
+ /// Index of column A
+ /// Index of column B
+ /// scalar cos value
+ /// scalar sin value
+ private static void Csrot(Matrix a, int rowCount, int columnA, int columnB, double c, double s)
+ {
+ for (var i = 0; i < rowCount; i++)
+ {
+ var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB));
+ var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA));
+ a.At(i, columnB, tmp);
+ a.At(i, columnA, z);
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, AX = B, with A SVD factorized.
+ ///
+ /// The right hand side , B.
+ /// The left hand side , X.
+ public override void Solve(Matrix input, Matrix result)
+ {
+ // Check for proper arguments.
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (!ComputeVectors)
+ {
+ throw new InvalidOperationException(Resources.SingularVectorsNotComputed);
+ }
+
+ // The solution X should have the same number of columns as B
+ if (input.ColumnCount != result.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
+ if (MatrixU.RowCount != input.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+ }
+
+ // The solution X row dimension is equal to the column dimension of A
+ if (MatrixVT.ColumnCount != result.RowCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+ }
+
+ var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount);
+ var bn = input.ColumnCount;
+
+ var tmp = new Complex[MatrixVT.ColumnCount];
+
+ for (var k = 0; k < bn; k++)
+ {
+ for (var j = 0; j < MatrixVT.ColumnCount; j++)
+ {
+ var value = Complex.Zero;
+ if (j < mn)
+ {
+ for (var i = 0; i < MatrixU.RowCount; i++)
+ {
+ value += MatrixU.At(i, j).Conjugate() * input.At(i, k);
+ }
+
+ value /= VectorS[j];
+ }
+
+ tmp[j] = value;
+ }
+
+ for (var j = 0; j < MatrixVT.ColumnCount; j++)
+ {
+ var value = Complex.Zero;
+ for (var i = 0; i < MatrixVT.ColumnCount; i++)
+ {
+ value += MatrixVT.At(i, j).Conjugate() * tmp[i];
+ }
+
+ result[j, k] = value;
+ }
+ }
+ }
+
+ ///
+ /// Solves a system of linear equations, Ax = b, with A SVD factorized.
+ ///
+ /// The right hand side vector, b.
+ /// The left hand side , x.
+ public override void Solve(Vector input, Vector result)
+ {
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (!ComputeVectors)
+ {
+ throw new InvalidOperationException(Resources.SingularVectorsNotComputed);
+ }
+
+ // Ax=b where A is an m x n matrix
+ // Check that b is a column vector with m entries
+ if (MatrixU.RowCount != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Check that x is a column vector with n entries
+ if (MatrixVT.ColumnCount != result.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount);
+ var tmp = new Complex[MatrixVT.ColumnCount];
+ for (var j = 0; j < MatrixVT.ColumnCount; j++)
+ {
+ var value = Complex.Zero;
+ if (j < mn)
+ {
+ for (var i = 0; i < MatrixU.RowCount; i++)
+ {
+ value += MatrixU.At(i, j).Conjugate() * input[i];
+ }
+
+ value /= VectorS[j];
+ }
+
+ tmp[j] = value;
+ }
+
+ for (var j = 0; j < MatrixVT.ColumnCount; j++)
+ {
+ var value = Complex.Zero;
+ for (var i = 0; i < MatrixVT.ColumnCount; i++)
+ {
+ value += MatrixVT.At(i, j).Conjugate() * tmp[i];
+ }
+
+ result[j] = value;
+ }
+ }
+
+ #region Simple arithmetic of type T
+
+ ///
+ /// Multiply two values T*T
+ ///
+ /// Left operand value
+ /// Right operand value
+ /// Result of multiplication
+ protected sealed override Complex MultiplyT(Complex val1, Complex val2)
+ {
+ return val1 * val2;
+ }
+
+ ///
+ /// Returns the absolute value of a specified number.
+ ///
+ /// A number whose absolute is to be found
+ /// Absolute value
+ protected sealed override double AbsoluteT(Complex val1)
+ {
+ return val1.Magnitude;
+ }
+
+ ///
+ /// Get value of type T equal to one
+ ///
+ /// One value
+ protected sealed override Complex OneValueT
+ {
+ get { return Complex.One; }
+ }
+
+ #endregion
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs
new file mode 100644
index 00000000..356eefca
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs
@@ -0,0 +1,521 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Solvers;
+ using Generic.Solvers.Preconditioners;
+ using Generic.Solvers.Status;
+ using Preconditioners;
+ using Properties;
+
+ ///
+ /// A Bi-Conjugate Gradient stabilized iterative matrix solver.
+ ///
+ ///
+ ///
+ /// The Bi-Conjugate Gradient Stabilized (BiCGStab) solver is an 'improvement'
+ /// of the standard Conjugate Gradient (CG) solver. Unlike the CG solver the
+ /// BiCGStab can be used on non-symmetric matrices.
+ /// Note that much of the success of the solver depends on the selection of the
+ /// proper preconditioner.
+ ///
+ ///
+ /// The Bi-CGSTAB algorithm was taken from:
+ /// Templates for the solution of linear systems: Building blocks
+ /// for iterative methods
+ ///
+ /// Richard Barrett, Michael Berry, Tony F. Chan, James Demmel,
+ /// June M. Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo,
+ /// Charles Romine and Henk van der Vorst
+ ///
+ /// Url: http://www.netlib.org/templates/Templates.html
+ ///
+ /// Algorithm is described in Chapter 2, section 2.3.8, page 27
+ ///
+ ///
+ /// The example code below provides an indication of the possible use of the
+ /// solver.
+ ///
+ ///
+ public sealed class BiCgStab : IIterativeSolver
+ {
+ ///
+ /// The status used if there is no status, i.e. the solver hasn't run yet and there is no
+ /// iterator.
+ ///
+ private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined();
+
+ ///
+ /// The preconditioner that will be used. Can be set to , in which case the default
+ /// pre-conditioner will be used.
+ ///
+ private IPreConditioner _preconditioner;
+
+ ///
+ /// The iterative process controller.
+ ///
+ private IIterator _iterator;
+
+ ///
+ /// Indicates if the user has stopped the solver.
+ ///
+ private bool _hasBeenStopped;
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ /// When using this constructor the solver will use the with
+ /// the standard settings and a default preconditioner.
+ ///
+ public BiCgStab() : this(null, null)
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ ///
+ /// When using this constructor the solver will use a default preconditioner.
+ ///
+ ///
+ /// The main advantages of using a user defined are:
+ ///
+ /// - It is possible to set the desired convergence limits.
+ /// -
+ /// It is possible to check the reason for which the solver finished
+ /// the iterative procedure by calling the property.
+ ///
+ ///
+ ///
+ ///
+ /// The that will be used to monitor the iterative process.
+ public BiCgStab(IIterator iterator) : this(null, iterator)
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ /// When using this constructor the solver will use the with
+ /// the standard settings.
+ ///
+ /// The that will be used to precondition the matrix equation.
+ public BiCgStab(IPreConditioner preconditioner) : this(preconditioner, null)
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ ///
+ /// The main advantages of using a user defined are:
+ ///
+ /// - It is possible to set the desired convergence limits.
+ /// -
+ /// It is possible to check the reason for which the solver finished
+ /// the iterative procedure by calling the property.
+ ///
+ ///
+ ///
+ ///
+ /// The that will be used to precondition the matrix equation.
+ /// The that will be used to monitor the iterative process.
+ public BiCgStab(IPreConditioner preconditioner, IIterator iterator)
+ {
+ _iterator = iterator;
+ _preconditioner = preconditioner;
+ }
+
+ ///
+ /// Sets the that will be used to precondition the iterative process.
+ ///
+ /// The preconditioner.
+ public void SetPreconditioner(IPreConditioner preconditioner)
+ {
+ _preconditioner = preconditioner;
+ }
+
+ ///
+ /// Sets the that will be used to track the iterative process.
+ ///
+ /// The iterator.
+ public void SetIterator(IIterator iterator)
+ {
+ _iterator = iterator;
+ }
+
+ ///
+ /// Gets the status of the iteration once the calculation is finished.
+ ///
+ public ICalculationStatus IterationResult
+ {
+ get
+ {
+ return (_iterator != null) ? _iterator.Status : DefaultStatus;
+ }
+ }
+
+ ///
+ /// Stops the solve process.
+ ///
+ ///
+ /// Note that it may take an indetermined amount of time for the solver to actually stop the process.
+ ///
+ public void StopSolve()
+ {
+ _hasBeenStopped = true;
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient , A.
+ /// The solution , b.
+ /// The result , x.
+ public Vector Solve(Matrix matrix, Vector vector)
+ {
+ if (vector == null)
+ {
+ throw new ArgumentNullException();
+ }
+
+ Vector result = new DenseVector(matrix.RowCount);
+ Solve(matrix, vector, result);
+ return result;
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient , A.
+ /// The solution , b.
+ /// The result , x.
+ public void Solve(Matrix matrix, Vector input, Vector result)
+ {
+ // If we were stopped before, we are no longer
+ // We're doing this at the start of the method to ensure
+ // that we can use these fields immediately.
+ _hasBeenStopped = false;
+
+ // Parameters checks
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
+ }
+
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.Count != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Initialize the solver fields
+ // Set the convergence monitor
+ if (_iterator == null)
+ {
+ _iterator = Iterator.CreateDefault();
+ }
+
+ if (_preconditioner == null)
+ {
+ _preconditioner = new UnitPreconditioner();
+ }
+
+ _preconditioner.Initialize(matrix);
+
+ // Compute r_0 = b - Ax_0 for some initial guess x_0
+ // In this case we take x_0 = vector
+ // This is basically a SAXPY so it could be made a lot faster
+ Vector residuals = new DenseVector(matrix.RowCount);
+ CalculateTrueResidual(matrix, residuals, result, input);
+
+ // Choose r~ (for example, r~ = r_0)
+ var tempResiduals = residuals.Clone();
+ var temp = result.Clone();
+
+ // create five temporary vectors needed to hold temporary
+ // coefficients. All vectors are mangled in each iteration.
+ // These are defined here to prevent stressing the garbage collector
+ Vector vecP = new DenseVector(residuals.Count);
+ Vector vecPdash = new DenseVector(residuals.Count);
+ Vector nu = new DenseVector(residuals.Count);
+ Vector vecS = new DenseVector(residuals.Count);
+ Vector vecSdash = new DenseVector(residuals.Count);
+ Vector t = new DenseVector(residuals.Count);
+
+ // create some temporary double variables that are needed
+ // to hold values in between iterations
+ Complex currentRho = 0;
+ Complex alpha = 0;
+ Complex omega = 0;
+
+ var iterationNumber = 0;
+ while (ShouldContinue(iterationNumber, result, input, residuals))
+ {
+ // rho_(i-1) = r~^T r_(i-1) // dotproduct r~ and r_(i-1)
+ var oldRho = currentRho;
+ currentRho = tempResiduals.DotProduct(residuals);
+
+ // if (rho_(i-1) == 0) // METHOD FAILS
+ // If rho is only 1 ULP from zero then we fail.
+ if (currentRho.Real.AlmostEqual(0, 1) && currentRho.Imaginary.AlmostEqual(0, 1))
+ {
+ // Rho-type breakdown
+ throw new Exception("Iterative solver experience a numerical break down");
+ }
+
+ if (iterationNumber != 0)
+ {
+ // beta_(i-1) = (rho_(i-1)/rho_(i-2))(alpha_(i-1)/omega(i-1))
+ var beta = (currentRho / oldRho) * (alpha / omega);
+
+ // p_i = r_(i-1) + beta_(i-1)(p_(i-1) - omega_(i-1) * nu_(i-1))
+ vecP.Add(nu.Multiply(-omega), vecP);
+
+ vecP.Multiply(beta, vecP);
+ vecP.Add(residuals, vecP);
+ }
+ else
+ {
+ // p_i = r_(i-1)
+ residuals.CopyTo(vecP);
+ }
+
+ // SOLVE Mp~ = p_i // M = preconditioner
+ _preconditioner.Approximate(vecP, vecPdash);
+
+ // nu_i = Ap~
+ matrix.Multiply(vecPdash, nu);
+
+ // alpha_i = rho_(i-1)/ (r~^T nu_i) = rho / dotproduct(r~ and nu_i)
+ alpha = currentRho * 1 / tempResiduals.DotProduct(nu);
+
+ // s = r_(i-1) - alpha_i nu_i
+ residuals.Add(nu.Multiply(-alpha), vecS);
+
+ // Check if we're converged. If so then stop. Otherwise continue;
+ // Calculate the temporary result.
+ // Be careful not to change any of the temp vectors, except for
+ // temp. Others will be used in the calculation later on.
+ // x_i = x_(i-1) + alpha_i * p^_i + s^_i
+ vecPdash.Multiply(alpha, temp);
+ temp.Add(vecSdash, temp);
+ temp.Add(result, temp);
+
+ // Check convergence and stop if we are converged.
+ if (!ShouldContinue(iterationNumber, temp, input, vecS))
+ {
+ temp.CopyTo(result);
+
+ // Calculate the true residual
+ CalculateTrueResidual(matrix, residuals, result, input);
+
+ // Now recheck the convergence
+ if (!ShouldContinue(iterationNumber, result, input, residuals))
+ {
+ // We're all good now.
+ return;
+ }
+
+ // Continue the calculation
+ iterationNumber++;
+ continue;
+ }
+
+ // SOLVE Ms~ = s
+ _preconditioner.Approximate(vecS, vecSdash);
+
+ // temp = As~
+ matrix.Multiply(vecSdash, t);
+
+ // omega_i = temp^T s / temp^T temp
+ omega = t.DotProduct(vecS) / t.DotProduct(t);
+
+ // x_i = x_(i-1) + alpha_i p^ + omega_i s^
+ result.Add(vecSdash.Multiply(omega), result);
+ result.Add(vecPdash.Multiply(alpha), result);
+
+ t.Multiply(-omega, residuals);
+ residuals.Add(vecS, residuals);
+
+ // for continuation it is necessary that omega_i != 0.0
+ // If omega is only 1 ULP from zero then we fail.
+ if (omega.Real.AlmostEqual(0, 1) && omega.Imaginary.AlmostEqual(0, 1))
+ {
+ // Omega-type breakdown
+ throw new Exception("Iterative solver experience a numerical break down");
+ }
+
+ if (!ShouldContinue(iterationNumber, result, input, residuals))
+ {
+ // Recalculate the residuals and go round again. This is done to ensure that
+ // we have the proper residuals.
+ // The residual calculation based on omega_i * s can be off by a factor 10. So here
+ // we calculate the real residual (which can be expensive) but we only do it if we're
+ // sufficiently close to the finish.
+ CalculateTrueResidual(matrix, residuals, result, input);
+ }
+
+ iterationNumber++;
+ }
+ }
+
+ ///
+ /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax
+ ///
+ /// Instance of the A.
+ /// Residual values in .
+ /// Instance of the x.
+ /// Instance of the b.
+ private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+ {
+ // -Ax = residual
+ matrix.Multiply(x, residual);
+
+ // Do not use residual = residual.Negate() because it creates another object
+ residual.Multiply(-1, residual);
+
+ // residual + b
+ residual.Add(b, residual);
+ }
+
+ ///
+ /// Determine if calculation should continue
+ ///
+ /// Number of iterations passed
+ /// Result .
+ /// Source .
+ /// Residual .
+ /// true if continue, otherwise false
+ private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+ {
+ if (_hasBeenStopped)
+ {
+ _iterator.IterationCancelled();
+ return true;
+ }
+
+ _iterator.DetermineStatus(iterationNumber, result, source, residuals);
+ var status = _iterator.Status;
+
+ // We stop if either:
+ // - the user has stopped the calculation
+ // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.)
+ return (!status.TerminatesCalculation) && (!_hasBeenStopped);
+ }
+
+ ///
+ /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
+ /// solution matrix and X is the unknown matrix.
+ ///
+ /// The coefficient , A.
+ /// The solution , B.
+ /// The result , X.
+ public Matrix Solve(Matrix matrix, Matrix input)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+ Solve(matrix, input, result);
+ return result;
+ }
+
+ ///
+ /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
+ /// solution matrix and X is the unknown matrix.
+ ///
+ /// The coefficient , A.
+ /// The solution , B.
+ /// The result , X
+ public void Solve(Matrix matrix, Matrix input, Matrix result)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ for (var column = 0; column < input.ColumnCount; column++)
+ {
+ var solution = Solve(matrix, input.Column(column));
+ foreach (var element in solution.GetIndexedEnumerator())
+ {
+ result.At(element.Key, column, element.Value);
+ }
+ }
+ }
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/CompositeSolver.cs
new file mode 100644
index 00000000..0f4fca20
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/CompositeSolver.cs
@@ -0,0 +1,630 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
+{
+ using System;
+ using System.Collections.Generic;
+ using System.IO;
+ using System.Linq;
+ using System.Numerics;
+ using System.Reflection;
+ using Generic;
+ using Generic.Solvers;
+ using Generic.Solvers.Status;
+ using Properties;
+
+ ///
+ /// A composite matrix solver. The actual solver is made by a sequence of
+ /// matrix solvers.
+ ///
+ ///
+ ///
+ /// Solver based on:
+ /// Faster PDE-based simulations using robust composite linear solvers
+ /// S. Bhowmicka, P. Raghavan a,*, L. McInnes b, B. Norris
+ /// Future Generation Computer Systems, Vol 20, 2004, pp 373–387
+ ///
+ ///
+ /// Note that if an iterator is passed to this solver it will be used for all the sub-solvers.
+ ///
+ ///
+ public sealed class CompositeSolver : IIterativeSolver
+ {
+ #region Internal class - DoubleComparer
+ ///
+ /// An IComparer used to compare double precision floating points.
+ ///
+ /// NOTE: The instance of this class is used only in . If C# suppports interface inheritence
+ /// NOTE: and methods in anonymous types, then this class should be deleted and anonymous type implemented with IComaprer support
+ /// NOTE: in constructor
+ public sealed class DoubleComparer : IComparer
+ {
+ ///
+ /// Compares two double values based on the selected comparison method.
+ ///
+ /// The first double to compare.
+ /// The second double to compare.
+ ///
+ /// A 32-bit signed integer that indicates the relative order of the objects being compared. The return
+ /// value has the following meanings:
+ /// Value Meaning Less than zero This object is less than the other parameter.
+ /// Zero This object is equal to other.
+ /// Greater than zero This object is greater than other.
+ ///
+ public int Compare(double x, double y)
+ {
+ return x.CompareTo(y, 1);
+ }
+ }
+ #endregion
+
+ ///
+ /// The default status used if the solver is not running.
+ ///
+ private static readonly ICalculationStatus NonRunningStatus = new CalculationIndetermined();
+
+ ///
+ /// The default status used if the solver is running.
+ ///
+ private static readonly ICalculationStatus RunningStatus = new CalculationRunning();
+
+#if SILVERLIGHT
+ private static readonly Dictionary>> SolverSetups = new Dictionary>>();
+#else
+ ///
+ /// The collection of iterative solver setups. Stored based on the
+ /// ratio between the relative speed and relative accuracy.
+ ///
+ private static readonly SortedList>> SolverSetups = new SortedList>>(new DoubleComparer());
+#endif
+
+ #region Solver information loading methods
+
+ ///
+ /// Loads all the available objects from the MathNet.Numerics assembly.
+ ///
+ public static void LoadSolverInformation()
+ {
+ LoadSolverInformation(new Type[0]);
+ }
+
+ ///
+ /// Loads the available objects from the MathNet.Numerics assembly.
+ ///
+ /// The types that should not be loaded.
+ public static void LoadSolverInformation(Type[] typesToExclude)
+ {
+ LoadSolverInformationFromAssembly(Assembly.GetExecutingAssembly(), typesToExclude);
+ }
+
+ ///
+ /// Loads the available objects from the assembly specified by the file location.
+ ///
+ /// The fully qualified path to the assembly.
+ public static void LoadSolverInformationFromAssembly(string assemblyLocation)
+ {
+ LoadSolverInformationFromAssembly(assemblyLocation, new Type[0]);
+ }
+
+ ///
+ /// Loads the available objects from the assembly specified by the file location.
+ ///
+ /// The fully qualified path to the assembly.
+ /// The types that should not be loaded.
+ public static void LoadSolverInformationFromAssembly(string assemblyLocation, params Type[] typesToExclude)
+ {
+ if (assemblyLocation == null)
+ {
+ throw new ArgumentNullException("assemblyLocation");
+ }
+
+ if (assemblyLocation.Length == 0)
+ {
+ throw new ArgumentException();
+ }
+
+ if (!File.Exists(assemblyLocation))
+ {
+ throw new FileNotFoundException();
+ }
+
+ // Get the assembly name
+ var assemblyFileName = Path.GetFileNameWithoutExtension(assemblyLocation);
+
+ // Now load the assembly with an AssemblyName
+ var assemblyName = new AssemblyName(assemblyFileName);
+ var assembly = Assembly.Load(assemblyName);
+
+ // Can't get this because we checked that the file exists.
+ // FileNotFoundException --> Can't get this because we checked that the file exists.
+ // FileLoadException
+ // BadImageFormatException
+ // Now we can load the solver information.
+ LoadSolverInformationFromAssembly(assembly, typesToExclude);
+ }
+
+ ///
+ /// Loads the available objects from the assembly specified by the assembly name.
+ ///
+ /// The of the assembly that should be searched for setup objects.
+ public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName)
+ {
+ LoadSolverInformationFromAssembly(assemblyName, new Type[0]);
+ }
+
+ ///
+ /// Loads the available objects from the assembly specified by the assembly name.
+ ///
+ /// The of the assembly that should be searched for setup objects.
+ /// The types that should not be loaded.
+ public static void LoadSolverInformationFromAssembly(AssemblyName assemblyName, params Type[] typesToExclude)
+ {
+ if (assemblyName == null)
+ {
+ throw new ArgumentNullException("assemblyName");
+ }
+
+ var assembly = Assembly.Load(assemblyName);
+
+ // May throw:
+ // ArgumentNullException --> Can't get this because we checked it.
+ // FileNotFoundException
+ // FileLoadException
+ // BadImageFormatException
+ // Now we can load the solver information.
+ LoadSolverInformationFromAssembly(assembly, typesToExclude);
+ }
+
+ ///
+ /// Loads the available objects from the assembly specified by the type.
+ ///
+ /// The type in the assembly which should be searched for setup objects.
+ public static void LoadSolverInformationFromAssembly(Type typeInAssembly)
+ {
+ LoadSolverInformationFromAssembly(typeInAssembly, new Type[0]);
+ }
+
+ ///
+ /// Loads the available objects from the assembly specified by the type.
+ ///
+ /// The type in the assembly which should be searched for setup objects.
+ /// The types that should not be loaded.
+ public static void LoadSolverInformationFromAssembly(Type typeInAssembly, params Type[] typesToExclude)
+ {
+ if (typeInAssembly == null)
+ {
+ throw new ArgumentNullException("typeInAssembly");
+ }
+
+ LoadSolverInformationFromAssembly(typeInAssembly.Assembly, typesToExclude);
+ }
+
+ ///
+ /// Loads the available objects from the specified assembly.
+ ///
+ /// The assembly which will be searched for setup objects.
+ public static void LoadSolverInformationFromAssembly(Assembly assembly)
+ {
+ LoadSolverInformationFromAssembly(assembly, new Type[0]);
+ }
+
+ ///
+ /// Loads the available objects from the specified assembly.
+ ///
+ /// The assembly which will be searched for setup objects.
+ /// The types that should not be loaded.
+ public static void LoadSolverInformationFromAssembly(Assembly assembly, params Type[] typesToExclude)
+ {
+ if (assembly == null)
+ {
+ throw new ArgumentNullException("Assembly");
+ }
+
+ if (typesToExclude == null)
+ {
+ throw new ArgumentNullException("typesToExclude");
+ }
+
+ var excludedTypes = new List(typesToExclude);
+
+ // Load all the types in the assembly
+ // Find all the types that implement IIterativeSolverSetup
+ // Create an object of each of these types
+ // Get the type of the iterative solver that will be instantiated by the setup object
+ // Check if it's on the excluding list, if so throw the setup object away otherwise keep it.
+ var interfaceTypes = new List();
+ foreach (var type in assembly.GetTypes().Where(type => (!type.IsAbstract && !type.IsEnum && !type.IsInterface && type.IsVisible)))
+ {
+ interfaceTypes.AddRange(type.GetInterfaces());
+ if (!interfaceTypes.Any(match => typeof(IIterativeSolverSetup).IsAssignableFrom(match)))
+ {
+ continue;
+ }
+
+ // See if we actually want this type of iterative solver
+ IIterativeSolverSetup setup;
+ try
+ {
+ // If something goes wrong we just ignore it and move on with the next type.
+ // There should probably be a log somewhere indicating that something went wrong?
+ setup = (IIterativeSolverSetup)Activator.CreateInstance(type);
+ }
+ catch (ArgumentException)
+ {
+ continue;
+ }
+ catch (NotSupportedException)
+ {
+ continue;
+ }
+ catch (TargetInvocationException)
+ {
+ continue;
+ }
+ catch (MethodAccessException)
+ {
+ continue;
+ }
+ catch (MissingMethodException)
+ {
+ continue;
+ }
+ catch (MemberAccessException)
+ {
+ continue;
+ }
+ catch (TypeLoadException)
+ {
+ continue;
+ }
+
+ if (excludedTypes.Any(match => match.IsAssignableFrom(setup.SolverType) ||
+ match.IsAssignableFrom(setup.PreconditionerType)))
+ {
+ continue;
+ }
+
+ // Ok we want the solver, so store the object
+ var ratio = setup.SolutionSpeed / setup.Reliability;
+ if (!SolverSetups.ContainsKey(ratio))
+ {
+ SolverSetups.Add(ratio, new List>());
+ }
+
+ var list = SolverSetups[ratio];
+ list.Add(setup);
+ }
+ }
+
+ #endregion
+
+ ///
+ /// The collection of solvers that will be used to
+ ///
+ private readonly List> _solvers = new List>();
+
+ ///
+ /// The status of the calculation.
+ ///
+ private ICalculationStatus _status = NonRunningStatus;
+
+ ///
+ /// The iterator that is used to control the iteration process.
+ ///
+ private IIterator _iterator;
+
+ ///
+ /// A flag indicating if the solver has been stopped or not.
+ ///
+ private bool _hasBeenStopped;
+
+ ///
+ /// The solver that is currently running. Reference is used to be able to stop the
+ /// solver if the user cancels the solve process.
+ ///
+ private IIterativeSolver _currentSolver;
+
+ ///
+ /// Initializes a new instance of the class with the default iterator.
+ ///
+ public CompositeSolver() : this(null)
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class with the specified iterator.
+ ///
+ /// The iterator that will be used to control the iteration process.
+ public CompositeSolver(IIterator iterator)
+ {
+ _iterator = iterator;
+ }
+
+ ///
+ /// Sets the IIterator that will be used to track the iterative process.
+ ///
+ /// The iterator.
+ public void SetIterator(IIterator iterator)
+ {
+ _iterator = iterator;
+ }
+
+ ///
+ /// Gets the status of the iteration once the calculation is finished.
+ ///
+ public ICalculationStatus IterationResult
+ {
+ get
+ {
+ return _status;
+ }
+ }
+
+ ///
+ /// Stops the solve process.
+ ///
+ ///
+ /// Note that it may take an indetermined amount of time for the solver to actually stop the process.
+ ///
+ public void StopSolve()
+ {
+ _hasBeenStopped = true;
+ if (_currentSolver != null)
+ {
+ _currentSolver.StopSolve();
+ }
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution vector, b.
+ /// The result vector, x.
+ public Vector Solve(Matrix matrix, Vector vector)
+ {
+ if (vector == null)
+ {
+ throw new ArgumentNullException();
+ }
+
+ Vector result = new DenseVector(matrix.RowCount);
+ Solve(matrix, vector, result);
+ return result;
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution vector, b
+ /// The result vector, x
+ public void Solve(Matrix matrix, Vector input, Vector result)
+ {
+ // If we were stopped before, we are no longer
+ // We're doing this at the start of the method to ensure
+ // that we can use these fields immediately.
+ _hasBeenStopped = false;
+ _currentSolver = null;
+
+ // Error checks
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
+ }
+
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.Count != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Initialize the solver fields
+ // Set the convergence monitor
+ if (_iterator == null)
+ {
+ _iterator = Iterator.CreateDefault();
+ }
+
+ // Load the solvers into our own internal data structure
+ // Once we have solvers we can always reuse them.
+ if (_solvers.Count == 0)
+ {
+ LoadSolvers();
+ }
+
+ // Create a copy of the solution and result vectors so we can use them
+ // later on
+ var internalInput = input.Clone();
+ var internalResult = result.Clone();
+
+ foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped))
+ {
+ // Store a reference to the solver so we can stop it.
+ _currentSolver = solver;
+
+ try
+ {
+ // Reset the iterator and pass it to the solver
+ _iterator.ResetToPrecalculationState();
+ solver.SetIterator(_iterator);
+
+ // Start the solver
+ solver.Solve(matrix, internalInput, internalResult);
+ }
+ catch (Exception)
+ {
+ // The solver broke down.
+ // Log a message about this
+ // Switch to the next preconditioner.
+ // Reset the solution vector to the previous solution
+ input.CopyTo(internalInput);
+ _status = RunningStatus;
+ continue;
+ }
+
+ // There was no fatal breakdown so check the status
+ if (_iterator.Status is CalculationConverged)
+ {
+ // We're done
+ break;
+ }
+
+ // We're not done
+ // Either:
+ // - calculation finished without convergence
+ if (_iterator.Status is CalculationStoppedWithoutConvergence)
+ {
+ // Copy the internal result to the result vector and
+ // continue with the calculation.
+ internalInput.CopyTo(input);
+ }
+ else
+ {
+ // - calculation failed --> restart with the original vector
+ // - calculation diverged --> restart with the original vector
+ // - Some unknown status occurred --> To be safe restart.
+ input.CopyTo(internalInput);
+ }
+ }
+
+ // Inside the loop we already copied the final results (if there are any)
+ // So no need to do that again.
+
+ // Clean up
+ // No longer need the current solver
+ _currentSolver = null;
+
+ // Set the final status
+ _status = _iterator.Status;
+ }
+
+ ///
+ /// Load solvers
+ ///
+ private void LoadSolvers()
+ {
+ if (SolverSetups.Count == 0)
+ {
+ throw new Exception("IIterativeSolverSetup objects not found");
+ }
+
+#if SILVERLIGHT
+ foreach (var setup in SolverSetups.OrderBy(solver => solver.Key, new DoubleComparer()).Select(pair => pair.Value).SelectMany(setups => setups))
+#else
+ foreach (var setup in SolverSetups.Select(pair => pair.Value).SelectMany(setups => setups))
+#endif
+ {
+ _solvers.Add(setup.CreateNew());
+ }
+ }
+
+ ///
+ /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
+ /// solution matrix and X is the unknown matrix.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution matrix, B.
+ /// The result matrix, X.
+ public Matrix Solve(Matrix matrix, Matrix input)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+ Solve(matrix, input, result);
+ return result;
+ }
+
+ ///
+ /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
+ /// solution matrix and X is the unknown matrix.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution matrix, B.
+ /// The result matrix, X
+ public void Solve(Matrix matrix, Matrix input, Matrix result)
+ {
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+ }
+
+ for (var column = 0; column < input.ColumnCount; column++)
+ {
+ var solution = Solve(matrix, input.Column(column));
+ foreach (var element in solution.GetIndexedEnumerator())
+ {
+ result.At(element.Key, column, element.Value);
+ }
+ }
+ }
+ }
+}
diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs
new file mode 100644
index 00000000..2b98d81e
--- /dev/null
+++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs
@@ -0,0 +1,624 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
+{
+ using System;
+ using System.Numerics;
+ using Generic;
+ using Generic.Solvers;
+ using Generic.Solvers.Preconditioners;
+ using Generic.Solvers.Status;
+ using Preconditioners;
+ using Properties;
+
+ ///
+ /// A Generalized Product Bi-Conjugate Gradient iterative matrix solver.
+ ///
+ ///
+ ///
+ /// The Generalized Product Bi-Conjugate Gradient (GPBiCG) solver is an
+ /// alternative version of the Bi-Conjugate Gradient stabilized (CG) solver.
+ /// Unlike the CG solver the GPBiCG solver can be used on
+ /// non-symmetric matrices.
+ /// Note that much of the success of the solver depends on the selection of the
+ /// proper preconditioner.
+ ///
+ ///
+ /// The GPBiCG algorithm was taken from:
+ /// GPBiCG(m,l): A hybrid of BiCGSTAB and GPBiCG methods with
+ /// efficiency and robustness
+ ///
+ /// S. Fujino
+ ///
+ /// Applied Numerical Mathematics, Volume 41, 2002, pp 107 - 117
+ ///
+ ///
+ ///
+ /// The example code below provides an indication of the possible use of the
+ /// solver.
+ ///
+ ///
+ public sealed class GpBiCg : IIterativeSolver
+ {
+ ///
+ /// The status used if there is no status, i.e. the solver hasn't run yet and there is no
+ /// iterator.
+ ///
+ private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined();
+
+ ///
+ /// The preconditioner that will be used. Can be set to null, in which case the default
+ /// pre-conditioner will be used.
+ ///
+ private IPreConditioner _preconditioner;
+
+ ///
+ /// The iterative process controller.
+ ///
+ private IIterator _iterator;
+
+ ///
+ /// Indicates the number of BiCGStab steps should be taken
+ /// before switching.
+ ///
+ private int _numberOfBiCgStabSteps = 1;
+
+ ///
+ /// Indicates the number of GPBiCG steps should be taken
+ /// before switching.
+ ///
+ private int _numberOfGpbiCgSteps = 4;
+
+ ///
+ /// Indicates if the user has stopped the solver.
+ ///
+ private bool _hasBeenStopped;
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ /// When using this constructor the solver will use the with
+ /// the standard settings and a default preconditioner.
+ ///
+ public GpBiCg() : this(null, null)
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ ///
+ /// When using this constructor the solver will use a default preconditioner.
+ ///
+ ///
+ /// The main advantages of using a user defined are:
+ ///
+ /// - It is possible to set the desired convergence limits.
+ /// -
+ /// It is possible to check the reason for which the solver finished
+ /// the iterative procedure by calling the property.
+ ///
+ ///
+ ///
+ ///
+ /// The that will be used to monitor the iterative process.
+ public GpBiCg(IIterator iterator) : this(null, iterator)
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ /// When using this constructor the solver will use the with
+ /// the standard settings.
+ ///
+ /// The that will be used to precondition the matrix equation.
+ public GpBiCg(IPreConditioner preconditioner) : this(preconditioner, null)
+ {
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ ///
+ ///
+ /// The main advantages of using a user defined are:
+ ///
+ /// - It is possible to set the desired convergence limits.
+ /// -
+ /// It is possible to check the reason for which the solver finished
+ /// the iterative procedure by calling the property.
+ ///
+ ///
+ ///
+ ///
+ /// The that will be used to precondition the matrix equation.
+ /// The that will be used to monitor the iterative process.
+ public GpBiCg(IPreConditioner preconditioner, IIterator iterator)
+ {
+ _iterator = iterator;
+ _preconditioner = preconditioner;
+ }
+
+ ///
+ /// Gets or sets the number of steps taken with the BiCgStab algorithm
+ /// before switching over to the GPBiCG algorithm.
+ ///
+ public int NumberOfBiCgStabSteps
+ {
+ get
+ {
+ return _numberOfBiCgStabSteps;
+ }
+
+ set
+ {
+ if (value < 0)
+ {
+ throw new ArgumentOutOfRangeException("value");
+ }
+
+ _numberOfBiCgStabSteps = value;
+ }
+ }
+
+ ///
+ /// Gets or sets the number of steps taken with the GPBiCG algorithm
+ /// before switching over to the BiCgStab algorithm.
+ ///
+ public int NumberOfGpBiCgSteps
+ {
+ get
+ {
+ return _numberOfGpbiCgSteps;
+ }
+
+ set
+ {
+ if (value < 0)
+ {
+ throw new ArgumentOutOfRangeException("value");
+ }
+
+ _numberOfGpbiCgSteps = value;
+ }
+ }
+
+ ///
+ /// Sets the that will be used to precondition the iterative process.
+ ///
+ /// The preconditioner.
+ public void SetPreconditioner(IPreConditioner preconditioner)
+ {
+ _preconditioner = preconditioner;
+ }
+
+ ///
+ /// Sets the that will be used to track the iterative process.
+ ///
+ /// The iterator.
+ public void SetIterator(IIterator iterator)
+ {
+ _iterator = iterator;
+ }
+
+ ///
+ /// Gets the status of the iteration once the calculation is finished.
+ ///
+ public ICalculationStatus IterationResult
+ {
+ get
+ {
+ return (_iterator != null) ? _iterator.Status : DefaultStatus;
+ }
+ }
+
+ ///
+ /// Stops the solve process.
+ ///
+ ///
+ /// Note that it may take an indetermined amount of time for the solver to actually
+ /// stop the process.
+ ///
+ public void StopSolve()
+ {
+ _hasBeenStopped = true;
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution vector, b.
+ /// The result vector, x.
+ public Vector Solve(Matrix matrix, Vector vector)
+ {
+ if (vector == null)
+ {
+ throw new ArgumentNullException();
+ }
+
+ Vector result = new DenseVector(matrix.RowCount);
+ Solve(matrix, vector, result);
+ return result;
+ }
+
+ ///
+ /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
+ /// solution vector and x is the unknown vector.
+ ///
+ /// The coefficient matrix, A.
+ /// The solution vector, b
+ /// The result vector, x
+ public void Solve(Matrix matrix, Vector input, Vector result)
+ {
+ // If we were stopped before, we are no longer
+ // We're doing this at the start of the method to ensure
+ // that we can use these fields immediately.
+ _hasBeenStopped = false;
+
+ // Error checks
+ if (matrix == null)
+ {
+ throw new ArgumentNullException("matrix");
+ }
+
+ if (matrix.RowCount != matrix.ColumnCount)
+ {
+ throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
+ }
+
+ if (input == null)
+ {
+ throw new ArgumentNullException("input");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (result.Count != input.Count)
+ {
+ throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+ }
+
+ // Initialize the solver fields
+
+ // Set the convergence monitor
+ if (_iterator == null)
+ {
+ _iterator = Iterator.CreateDefault();
+ }
+
+ if (_preconditioner == null)
+ {
+ _preconditioner = new UnitPreconditioner();
+ }
+
+ _preconditioner.Initialize(matrix);
+
+ // x_0 is initial guess
+ // Take x_0 = 0
+ Vector