53 changed files with 9377 additions and 459 deletions
@ -0,0 +1,399 @@ |
|||||
|
// <copyright file="Poisson.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.Distributions |
||||
|
{ |
||||
|
using System; |
||||
|
using System.Collections.Generic; |
||||
|
using Properties; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Pseudo-random generation of poisson distributed deviates.
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// <para>Distribution is described at <a href="http://en.wikipedia.org/wiki/Poisson_distribution"> Wikipedia - Poisson distribution</a>.</para>
|
||||
|
/// <para>Knuth's method is used to generate Poisson distributed random variables.</para>
|
||||
|
/// <para>f(x) = exp(-λ)*λ^x/x!;</para>
|
||||
|
/// </remarks>
|
||||
|
public class Poisson : IDiscreteDistribution |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// The Poisson distribution parameter λ.
|
||||
|
/// </summary>
|
||||
|
private double _lambda; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// The distribution's random number generator.
|
||||
|
/// </summary>
|
||||
|
private Random _random; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets or sets the Poisson distribution parameter λ.
|
||||
|
/// </summary>
|
||||
|
public double Lambda |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return _lambda; |
||||
|
} |
||||
|
|
||||
|
set |
||||
|
{ |
||||
|
SetParameters(value); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="Poisson"/> class.
|
||||
|
/// </summary>
|
||||
|
/// <param name="lambda">The Poisson distribution parameter λ.</param>
|
||||
|
/// <exception cref="System.ArgumentOutOfRangeException">If <paramref name="lambda"/> is equal or less then 0.0.</exception>
|
||||
|
public Poisson(double lambda) |
||||
|
{ |
||||
|
SetParameters(lambda); |
||||
|
RandomSource = new Random(); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Sets the parameters of the distribution after checking their validity.
|
||||
|
/// </summary>
|
||||
|
/// <param name="lambda">The mean (λ) of the distribution.</param>
|
||||
|
/// <exception cref="System.ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
|
||||
|
private void SetParameters(double lambda) |
||||
|
{ |
||||
|
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) |
||||
|
{ |
||||
|
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); |
||||
|
} |
||||
|
|
||||
|
_lambda = lambda; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Checks whether the parameters of the distribution are valid.
|
||||
|
/// </summary>
|
||||
|
/// <param name="lambda">The mean (λ) of the distribution.</param>
|
||||
|
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
|
||||
|
private static bool IsValidParameterSet(double lambda) |
||||
|
{ |
||||
|
return lambda > 0.00; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Returns a <see cref="System.String"/> that represents this instance.
|
||||
|
/// </summary>
|
||||
|
/// <returns>
|
||||
|
/// A <see cref="System.String"/> that represents this instance.
|
||||
|
/// </returns>
|
||||
|
public override string ToString() |
||||
|
{ |
||||
|
return "Poisson(λ = " + _lambda + ")"; |
||||
|
} |
||||
|
|
||||
|
#region IDistribution Members
|
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets or sets the random number generator which is used to draw random samples.
|
||||
|
/// </summary>
|
||||
|
public Random RandomSource |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return _random; |
||||
|
} |
||||
|
|
||||
|
set |
||||
|
{ |
||||
|
if (value == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException(); |
||||
|
} |
||||
|
|
||||
|
_random = value; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the mean of the distribution.
|
||||
|
/// </summary>
|
||||
|
public double Mean |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return _lambda; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the variance of the distribution.
|
||||
|
/// </summary>
|
||||
|
public double Variance |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return _lambda; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the standard deviation of the distribution.
|
||||
|
/// </summary>
|
||||
|
public double StdDev |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return Math.Sqrt(_lambda); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the entropy of the distribution.
|
||||
|
/// </summary>
|
||||
|
/// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks>
|
||||
|
public double Entropy |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return (0.5 * Math.Log(2 * Constants.Pi * Constants.E * _lambda)) - (1.0 / (12.0 * _lambda)) - (1.0 / (24.0 * _lambda * _lambda)) - (19.0 / (360.0 * _lambda * _lambda * _lambda)); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the skewness of the distribution.
|
||||
|
/// </summary>
|
||||
|
public double Skewness |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return 1.0 / Math.Sqrt(_lambda); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the smallest element in the domain of the distributions which can be represented by an integer.
|
||||
|
/// </summary>
|
||||
|
public int Minimum |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return 0; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the largest element in the domain of the distributions which can be represented by an integer.
|
||||
|
/// </summary>
|
||||
|
public int Maximum |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return int.MaxValue; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Computes the cumulative distribution function of the Poisson distribution.
|
||||
|
/// </summary>
|
||||
|
/// <param name="x">The location at which to compute the cumulative density.</param>
|
||||
|
/// <returns>the cumulative density at <paramref name="x"/>.</returns>
|
||||
|
public double CumulativeDistribution(double x) |
||||
|
{ |
||||
|
return 1.0 - SpecialFunctions.GammaLowerRegularized(x + 1, _lambda); |
||||
|
} |
||||
|
|
||||
|
#endregion
|
||||
|
|
||||
|
#region IDiscreteDistribution Members
|
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the mode of the distribution.
|
||||
|
/// </summary>
|
||||
|
public int Mode |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return (int)Math.Floor(_lambda); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the median of the distribution.
|
||||
|
/// </summary>
|
||||
|
/// <remarks>Approximation, see Wikipedia <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a></remarks>
|
||||
|
public int Median |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return (int)Math.Floor(_lambda + (1.0 / 3.0) - (0.02 / _lambda)); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Computes values of the probability mass function.
|
||||
|
/// </summary>
|
||||
|
/// <param name="k">The location in the domain where we want to evaluate the probability mass function.</param>
|
||||
|
/// <returns>the probability mass at location <paramref name="k"/>.</returns>
|
||||
|
public double Probability(int k) |
||||
|
{ |
||||
|
return Math.Exp(-_lambda + (k * Math.Log(_lambda)) - SpecialFunctions.FactorialLn(k)); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Computes values of the log probability mass function.
|
||||
|
/// </summary>
|
||||
|
/// <param name="k">The location in the domain where we want to evaluate the log probability mass function.</param>
|
||||
|
/// <returns>the log probability mass at location <paramref name="k"/>.</returns>
|
||||
|
public double ProbabilityLn(int k) |
||||
|
{ |
||||
|
return -_lambda + (k * Math.Log(_lambda)) - SpecialFunctions.FactorialLn(k); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Samples a Poisson distributed random variable.
|
||||
|
/// </summary>
|
||||
|
/// <returns>A sample from the Poisson distribution.</returns>
|
||||
|
public int Sample() |
||||
|
{ |
||||
|
return DoSample(RandomSource, _lambda); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Samples an array of Poisson distributed random variables.
|
||||
|
/// </summary>
|
||||
|
/// <returns>a sequence of successes in N trials.</returns>
|
||||
|
public IEnumerable<int> Samples() |
||||
|
{ |
||||
|
while (true) |
||||
|
{ |
||||
|
yield return DoSample(RandomSource, _lambda); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
#endregion
|
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Samples a Poisson distributed random variable.
|
||||
|
/// </summary>
|
||||
|
/// <param name="rnd">The random number generator to use.</param>
|
||||
|
/// <param name="lambda">The Poisson distribution parameter λ.</param>
|
||||
|
/// <returns>A sample from the Poisson distribution.</returns>
|
||||
|
public static int Sample(Random rnd, double lambda) |
||||
|
{ |
||||
|
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) |
||||
|
{ |
||||
|
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); |
||||
|
} |
||||
|
|
||||
|
return DoSample(rnd, lambda); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Samples a sequence of Poisson distributed random variables.
|
||||
|
/// </summary>
|
||||
|
/// <param name="rnd">The random number generator to use.</param>
|
||||
|
/// <param name="lambda">The Poisson distribution parameter λ.</param>
|
||||
|
/// <returns>a sequence of samples from the distribution.</returns>
|
||||
|
public static IEnumerable<int> Samples(Random rnd, double lambda) |
||||
|
{ |
||||
|
if (Control.CheckDistributionParameters && !IsValidParameterSet(lambda)) |
||||
|
{ |
||||
|
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); |
||||
|
} |
||||
|
|
||||
|
while (true) |
||||
|
{ |
||||
|
yield return DoSample(rnd, lambda); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Generates one sample from the Poisson distribution.
|
||||
|
/// </summary>
|
||||
|
/// <param name="rnd">The random source to use.</param>
|
||||
|
/// <param name="lambda">The Poisson distribution parameter λ.</param>
|
||||
|
/// <returns>A random sample from the Poisson distribution.</returns>
|
||||
|
private static int DoSample(Random rnd, double lambda) |
||||
|
{ |
||||
|
return (lambda < 30.0) ? DoSampleShort(rnd, lambda) : DoSampleLarge(rnd, lambda); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Generates one sample from the Poisson distribution by Knuth's method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="rnd">The random source to use.</param>
|
||||
|
/// <param name="lambda">The Poisson distribution parameter λ.</param>
|
||||
|
/// <returns>A random sample from the Poisson distribution.</returns>
|
||||
|
private static int DoSampleShort(Random rnd, double lambda) |
||||
|
{ |
||||
|
var limit = Math.Exp(-lambda); |
||||
|
var count = 0; |
||||
|
for (var product = rnd.NextDouble(); product >= limit; product *= rnd.NextDouble()) |
||||
|
{ |
||||
|
count++; |
||||
|
} |
||||
|
|
||||
|
return count; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Generates one sample from the Poisson distribution by "Rejection method PA".
|
||||
|
/// </summary>
|
||||
|
/// <param name="rnd">The random source to use.</param>
|
||||
|
/// <param name="lambda">The Poisson distribution parameter λ.</param>
|
||||
|
/// <returns>A random sample from the Poisson distribution.</returns>
|
||||
|
/// <remarks>"Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson,
|
||||
|
/// Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
|
||||
|
/// The article is on pages 29-35. The algorithm given here is on page 32. </remarks>
|
||||
|
private static int DoSampleLarge(Random rnd, double lambda) |
||||
|
{ |
||||
|
var c = 0.767 - (3.36 / lambda); |
||||
|
var beta = Math.PI / Math.Sqrt(3.0 * lambda); |
||||
|
var alpha = beta * lambda; |
||||
|
var k = Math.Log(c) - lambda - Math.Log(beta); |
||||
|
|
||||
|
for (;;) |
||||
|
{ |
||||
|
var u = rnd.NextDouble(); |
||||
|
var x = (alpha - Math.Log((1.0 - u) / u)) / beta; |
||||
|
var n = (int)Math.Floor(x + 0.5); |
||||
|
if (n < 0) |
||||
|
{ |
||||
|
continue; |
||||
|
} |
||||
|
|
||||
|
var v = rnd.NextDouble(); |
||||
|
var y = alpha - (beta * x); |
||||
|
var temp = 1.0 + Math.Exp(y); |
||||
|
var lhs = y + Math.Log(v / (temp * temp)); |
||||
|
var rhs = k + (n * Math.Log(lambda)) - SpecialFunctions.FactorialLn(n); |
||||
|
if (lhs <= rhs) |
||||
|
{ |
||||
|
return n; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,964 @@ |
|||||
|
// <copyright file="DenseEvd.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using System.Numerics; |
||||
|
using Generic; |
||||
|
using Generic.Factorization; |
||||
|
using Properties; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Eigenvalues and eigenvectors of a complex matrix.
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// If A is hermitan, then A = V*D*V' where the eigenvalue matrix D is
|
||||
|
/// diagonal and the eigenvector matrix V is hermitan.
|
||||
|
/// I.e. A = V*D*V' and V*VH=I.
|
||||
|
/// If A is not symmetric, then the eigenvalue matrix D is block diagonal
|
||||
|
/// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
||||
|
/// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
|
||||
|
/// columns of V represent the eigenvectors in the sense that A*V = V*D,
|
||||
|
/// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly
|
||||
|
/// conditioned, or even singular, so the validity of the equation
|
||||
|
/// A = V*D*Inverse(V) depends upon V.cond().
|
||||
|
/// </remarks>
|
||||
|
public class DenseEvd : Evd<Complex> |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="DenseEvd"/> class. This object will compute the
|
||||
|
/// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
|
||||
|
/// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
|
||||
|
public DenseEvd(DenseMatrix matrix) |
||||
|
{ |
||||
|
if (matrix == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("matrix"); |
||||
|
} |
||||
|
|
||||
|
if (matrix.RowCount != matrix.ColumnCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
||||
|
} |
||||
|
|
||||
|
var order = matrix.RowCount; |
||||
|
|
||||
|
// Initialize matricies for eigenvalues and eigenvectors
|
||||
|
MatrixEv = DenseMatrix.Identity(order); |
||||
|
MatrixD = matrix.CreateMatrix(order, order); |
||||
|
VectorEv = new DenseVector(order); |
||||
|
|
||||
|
IsSymmetric = true; |
||||
|
|
||||
|
for (var i = 0; i < order & IsSymmetric; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < order & IsSymmetric; j++) |
||||
|
{ |
||||
|
IsSymmetric &= matrix[i, j] == matrix[j, i].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (IsSymmetric) |
||||
|
{ |
||||
|
var matrixCopy = matrix.ToArray(); |
||||
|
var tau = new Complex[order]; |
||||
|
var d = new double[order]; |
||||
|
var e = new double[order]; |
||||
|
|
||||
|
SymmetricTridiagonalize(matrixCopy, d, e, tau, order); |
||||
|
SymmetricDiagonalize(((DenseMatrix)MatrixEv).Data, d, e, order); |
||||
|
SymmetricUntridiagonalize(((DenseMatrix)MatrixEv).Data, matrixCopy, tau, order); |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
VectorEv[i] = new Complex(d[i], e[i]); |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
var matrixH = matrix.ToArray(); |
||||
|
NonsymmetricReduceToHessenberg(((DenseMatrix)MatrixEv).Data, matrixH, order); |
||||
|
NonsymmetricReduceHessenberToRealSchur(((DenseVector)VectorEv).Data, ((DenseMatrix)MatrixEv).Data, matrixH, order); |
||||
|
} |
||||
|
|
||||
|
MatrixD.SetDiagonal(VectorEv); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Reduces a complex hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrixA">Source matrix to reduce</param>
|
||||
|
/// <param name="d">Output: Arrays for internal storage of real parts of eigenvalues</param>
|
||||
|
/// <param name="e">Output: Arrays for internal storage of imaginary parts of eigenvalues</param>
|
||||
|
/// <param name="tau">Output: Arrays that contains further information about the transformations.</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures HTRIDI by
|
||||
|
/// Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for
|
||||
|
/// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void SymmetricTridiagonalize(Complex[,] matrixA, double[] d, double[] e, Complex[] tau, int order) |
||||
|
{ |
||||
|
double hh; |
||||
|
tau[order - 1] = Complex.One; |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
d[i] = matrixA[i, i].Real; |
||||
|
} |
||||
|
|
||||
|
// Householder reduction to tridiagonal form.
|
||||
|
for (var i = order - 1; i > 0; i--) |
||||
|
{ |
||||
|
// Scale to avoid under/overflow.
|
||||
|
var scale = 0.0; |
||||
|
var h = 0.0; |
||||
|
|
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
scale = scale + Math.Abs(matrixA[i, k].Real) + Math.Abs(matrixA[i, k].Imaginary); |
||||
|
} |
||||
|
|
||||
|
if (scale == 0.0) |
||||
|
{ |
||||
|
tau[i - 1] = Complex.One; |
||||
|
e[i] = 0.0; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
matrixA[i, k] /= scale; |
||||
|
h += matrixA[i, k].MagnitudeSquared(); |
||||
|
} |
||||
|
|
||||
|
Complex g = Math.Sqrt(h); |
||||
|
e[i] = scale * g.Real; |
||||
|
|
||||
|
Complex temp; |
||||
|
var f = matrixA[i, i - 1]; |
||||
|
if (f.Magnitude != 0) |
||||
|
{ |
||||
|
temp = -(matrixA[i, i - 1].Conjugate() * tau[i].Conjugate()) / f.Magnitude; |
||||
|
h += f.Magnitude * g.Real; |
||||
|
g = 1.0 + (g / f.Magnitude); |
||||
|
matrixA[i, i - 1] *= g; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
temp = -tau[i].Conjugate(); |
||||
|
matrixA[i, i - 1] = g; |
||||
|
} |
||||
|
|
||||
|
if ((f.Magnitude == 0) || (i != 1)) |
||||
|
{ |
||||
|
f = Complex.Zero; |
||||
|
for (var j = 0; j < i; j++) |
||||
|
{ |
||||
|
var tmp = Complex.Zero; |
||||
|
|
||||
|
// Form element of A*U.
|
||||
|
for (var k = 0; k <= j; k++) |
||||
|
{ |
||||
|
tmp += matrixA[j, k] * matrixA[i, k].Conjugate(); |
||||
|
} |
||||
|
|
||||
|
for (var k = j + 1; k <= i - 1; k++) |
||||
|
{ |
||||
|
tmp += matrixA[k, j].Conjugate() * matrixA[i, k].Conjugate(); |
||||
|
} |
||||
|
|
||||
|
// Form element of P
|
||||
|
tau[j] = tmp / h; |
||||
|
f += (tmp / h) * matrixA[i, j]; |
||||
|
} |
||||
|
|
||||
|
hh = f.Real / (h + h); |
||||
|
|
||||
|
// Form the reduced A.
|
||||
|
for (var j = 0; j < i; j++) |
||||
|
{ |
||||
|
f = matrixA[i, j].Conjugate(); |
||||
|
g = tau[j] - (hh * f); |
||||
|
tau[j] = g.Conjugate(); |
||||
|
|
||||
|
for (var k = 0; k <= j; k++) |
||||
|
{ |
||||
|
matrixA[j, k] -= (f * tau[k]) + (g * matrixA[i, k]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
matrixA[i, k] *= scale; |
||||
|
} |
||||
|
|
||||
|
tau[i - 1] = temp.Conjugate(); |
||||
|
} |
||||
|
|
||||
|
hh = d[i]; |
||||
|
d[i] = matrixA[i, i].Real; |
||||
|
matrixA[i, i] = new Complex(hh, scale * Math.Sqrt(h)); |
||||
|
} |
||||
|
|
||||
|
hh = d[0]; |
||||
|
d[0] = matrixA[0, 0].Real; |
||||
|
matrixA[0, 0] = hh; |
||||
|
e[0] = 0.0; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Symmetric tridiagonal QL algorithm.
|
||||
|
/// </summary>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="d">Arrays for internal storage of real parts of eigenvalues</param>
|
||||
|
/// <param name="e">Arrays for internal storage of imaginary parts of eigenvalues</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures tql2, by
|
||||
|
/// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
|
||||
|
/// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void SymmetricDiagonalize(Complex[] dataEv, double[] d, double[] e, int order) |
||||
|
{ |
||||
|
const int Maxiter = 1000; |
||||
|
|
||||
|
for (var i = 1; i < order; i++) |
||||
|
{ |
||||
|
e[i - 1] = e[i]; |
||||
|
} |
||||
|
|
||||
|
e[order - 1] = 0.0; |
||||
|
|
||||
|
var f = 0.0; |
||||
|
var tst1 = 0.0; |
||||
|
var eps = Precision.DoubleMachinePrecision; |
||||
|
for (var l = 0; l < order; l++) |
||||
|
{ |
||||
|
// Find small subdiagonal element
|
||||
|
tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l])); |
||||
|
var m = l; |
||||
|
while (m < order) |
||||
|
{ |
||||
|
if (Math.Abs(e[m]) <= eps * tst1) |
||||
|
{ |
||||
|
break; |
||||
|
} |
||||
|
|
||||
|
m++; |
||||
|
} |
||||
|
|
||||
|
// If m == l, d[l] is an eigenvalue,
|
||||
|
// otherwise, iterate.
|
||||
|
if (m > l) |
||||
|
{ |
||||
|
var iter = 0; |
||||
|
do |
||||
|
{ |
||||
|
iter = iter + 1; // (Could check iteration count here.)
|
||||
|
|
||||
|
// Compute implicit shift
|
||||
|
var g = d[l]; |
||||
|
var p = (d[l + 1] - g) / (2.0 * e[l]); |
||||
|
var r = SpecialFunctions.Hypotenuse(p, 1.0); |
||||
|
if (p < 0) |
||||
|
{ |
||||
|
r = -r; |
||||
|
} |
||||
|
|
||||
|
d[l] = e[l] / (p + r); |
||||
|
d[l + 1] = e[l] * (p + r); |
||||
|
|
||||
|
var dl1 = d[l + 1]; |
||||
|
var h = g - d[l]; |
||||
|
for (var i = l + 2; i < order; i++) |
||||
|
{ |
||||
|
d[i] -= h; |
||||
|
} |
||||
|
|
||||
|
f = f + h; |
||||
|
|
||||
|
// Implicit QL transformation.
|
||||
|
p = d[m]; |
||||
|
var c = 1.0; |
||||
|
var c2 = c; |
||||
|
var c3 = c; |
||||
|
var el1 = e[l + 1]; |
||||
|
var s = 0.0; |
||||
|
var s2 = 0.0; |
||||
|
for (var i = m - 1; i >= l; i--) |
||||
|
{ |
||||
|
c3 = c2; |
||||
|
c2 = c; |
||||
|
s2 = s; |
||||
|
g = c * e[i]; |
||||
|
h = c * p; |
||||
|
r = SpecialFunctions.Hypotenuse(p, e[i]); |
||||
|
e[i + 1] = s * r; |
||||
|
s = e[i] / r; |
||||
|
c = p / r; |
||||
|
p = (c * d[i]) - (s * g); |
||||
|
d[i + 1] = h + (s * ((c * g) + (s * d[i]))); |
||||
|
|
||||
|
// Accumulate transformation.
|
||||
|
for (var k = 0; k < order; k++) |
||||
|
{ |
||||
|
h = dataEv[((i + 1) * order) + k].Real; |
||||
|
dataEv[((i + 1) * order) + k] = (s * dataEv[(i * order) + k].Real) + (c * h); |
||||
|
dataEv[(i * order) + k] = (c * dataEv[(i * order) + k].Real) - (s * h); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
p = (-s) * s2 * c3 * el1 * e[l] / dl1; |
||||
|
e[l] = s * p; |
||||
|
d[l] = c * p; |
||||
|
|
||||
|
// Check for convergence. If too many iterations have been performed,
|
||||
|
// throw exception that Convergence Failed
|
||||
|
if (iter >= Maxiter) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ConvergenceFailed); |
||||
|
} |
||||
|
} |
||||
|
while (Math.Abs(e[l]) > eps * tst1); |
||||
|
} |
||||
|
|
||||
|
d[l] = d[l] + f; |
||||
|
e[l] = 0.0; |
||||
|
} |
||||
|
|
||||
|
// Sort eigenvalues and corresponding vectors.
|
||||
|
for (var i = 0; i < order - 1; i++) |
||||
|
{ |
||||
|
var k = i; |
||||
|
var p = d[i]; |
||||
|
for (var j = i + 1; j < order; j++) |
||||
|
{ |
||||
|
if (d[j] < p) |
||||
|
{ |
||||
|
k = j; |
||||
|
p = d[j]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (k != i) |
||||
|
{ |
||||
|
d[k] = d[i]; |
||||
|
d[i] = p; |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
p = dataEv[(i * order) + j].Real; |
||||
|
dataEv[(i * order) + j] = dataEv[(k * order) + j]; |
||||
|
dataEv[(k * order) + j] = p; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Determines eigenvectors by undoing the symmetric tridiagonalize transformation
|
||||
|
/// </summary>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="matrixA">Previously tridiagonalized matrix by <see cref="SymmetricTridiagonalize"/>.</param>
|
||||
|
/// <param name="tau">Contains further information about the transformations</param>
|
||||
|
/// <param name="order">Input matrix order</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures HTRIBK, by
|
||||
|
/// by Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for
|
||||
|
/// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void SymmetricUntridiagonalize(Complex[] dataEv, Complex[,] matrixA, Complex[] tau, int order) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
dataEv[(j * order) + i] = dataEv[(j * order) + i].Real * tau[i].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Recover and apply the Householder matrices.
|
||||
|
for (var i = 1; i < order; i++) |
||||
|
{ |
||||
|
var h = matrixA[i, i].Imaginary; |
||||
|
if (h != 0) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
var s = Complex.Zero; |
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
s += dataEv[(j * order) + k] * matrixA[i, k]; |
||||
|
} |
||||
|
|
||||
|
s = (s / h) / h; |
||||
|
|
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
dataEv[(j * order) + k] -= s * matrixA[i, k].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Nonsymmetric reduction to Hessenberg form.
|
||||
|
/// </summary>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="matrixH">Array for internal storage of nonsymmetric Hessenberg form.</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures orthes and ortran,
|
||||
|
/// by Martin and Wilkinson, Handbook for Auto. Comp.,
|
||||
|
/// Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutines in EISPACK.</remarks>
|
||||
|
private static void NonsymmetricReduceToHessenberg(Complex[] dataEv, Complex[,] matrixH, int order) |
||||
|
{ |
||||
|
var ort = new Complex[order]; |
||||
|
|
||||
|
for (var m = 1; m < order - 1; m++) |
||||
|
{ |
||||
|
// Scale column.
|
||||
|
var scale = 0.0; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
scale += Math.Abs(matrixH[i, m - 1].Real) + Math.Abs(matrixH[i, m - 1].Imaginary); |
||||
|
} |
||||
|
|
||||
|
if (scale != 0.0) |
||||
|
{ |
||||
|
// Compute Householder transformation.
|
||||
|
var h = 0.0; |
||||
|
for (var i = order - 1; i >= m; i--) |
||||
|
{ |
||||
|
ort[i] = matrixH[i, m - 1] / scale; |
||||
|
h += ort[i].MagnitudeSquared(); |
||||
|
} |
||||
|
|
||||
|
var g = Math.Sqrt(h); |
||||
|
if (ort[m].Magnitude != 0) |
||||
|
{ |
||||
|
h = h + (ort[m].Magnitude * g); |
||||
|
g /= ort[m].Magnitude; |
||||
|
ort[m] = (1.0 + g) * ort[m]; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
ort[m] = g; |
||||
|
matrixH[m, m - 1] = scale; |
||||
|
} |
||||
|
|
||||
|
// Apply Householder similarity transformation
|
||||
|
// H = (I-u*u'/h)*H*(I-u*u')/h)
|
||||
|
for (var j = m; j < order; j++) |
||||
|
{ |
||||
|
var f = Complex.Zero; |
||||
|
for (var i = order - 1; i >= m; i--) |
||||
|
{ |
||||
|
f += ort[i].Conjugate() * matrixH[i, j]; |
||||
|
} |
||||
|
|
||||
|
f = f / h; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
matrixH[i, j] -= f * ort[i]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
var f = Complex.Zero; |
||||
|
for (var j = order - 1; j >= m; j--) |
||||
|
{ |
||||
|
f += ort[j] * matrixH[i, j]; |
||||
|
} |
||||
|
|
||||
|
f = f / h; |
||||
|
for (var j = m; j < order; j++) |
||||
|
{ |
||||
|
matrixH[i, j] -= f * ort[j].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
ort[m] = scale * ort[m]; |
||||
|
matrixH[m, m - 1] *= -g; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Accumulate transformations (Algol's ortran).
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
dataEv[(j * order) + i] = i == j ? Complex.One : Complex.Zero; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var m = order - 2; m >= 1; m--) |
||||
|
{ |
||||
|
if (matrixH[m, m - 1] != Complex.Zero && ort[m] != Complex.Zero) |
||||
|
{ |
||||
|
var norm = (matrixH[m, m - 1].Real * ort[m].Real) + (matrixH[m, m - 1].Imaginary * ort[m].Imaginary); |
||||
|
|
||||
|
for (var i = m + 1; i < order; i++) |
||||
|
{ |
||||
|
ort[i] = matrixH[i, m - 1]; |
||||
|
} |
||||
|
|
||||
|
for (var j = m; j < order; j++) |
||||
|
{ |
||||
|
var g = Complex.Zero; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
g += ort[i].Conjugate() * dataEv[(j * order) + i]; |
||||
|
} |
||||
|
|
||||
|
// Double division avoids possible underflow
|
||||
|
g /= norm; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
dataEv[(j * order) + i] += g * ort[i]; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Create real subdiagonal elements.
|
||||
|
for (var i = 1; i < order; i++) |
||||
|
{ |
||||
|
if (matrixH[i, i - 1].Imaginary != 0.0) |
||||
|
{ |
||||
|
var y = matrixH[i, i - 1] / matrixH[i, i - 1].Magnitude; |
||||
|
matrixH[i, i - 1] = matrixH[i, i - 1].Magnitude; |
||||
|
for (var j = i; j < order; j++) |
||||
|
{ |
||||
|
matrixH[i, j] *= y.Conjugate(); |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j <= Math.Min(i + 1, order - 1); j++) |
||||
|
{ |
||||
|
matrixH[j, i] *= y; |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
dataEv[(i * order) + j] *= y; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Nonsymmetric reduction from Hessenberg to real Schur form.
|
||||
|
/// </summary>
|
||||
|
/// <param name="vectorV">Data array of the eigenvectors</param>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="matrixH">Array for internal storage of nonsymmetric Hessenberg form.</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedure hqr2,
|
||||
|
/// by Martin and Wilkinson, Handbook for Auto. Comp.,
|
||||
|
/// Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, Complex[] dataEv, Complex[,] matrixH, int order) |
||||
|
{ |
||||
|
// Initialize
|
||||
|
var n = order - 1; |
||||
|
var eps = Precision.DoubleMachinePrecision; |
||||
|
|
||||
|
double norm; |
||||
|
Complex s, x, y, z, exshift = Complex.Zero; |
||||
|
|
||||
|
// Outer loop over eigenvalue index
|
||||
|
var iter = 0; |
||||
|
while (n >= 0) |
||||
|
{ |
||||
|
// Look for single small sub-diagonal element
|
||||
|
var l = n; |
||||
|
while (l > 0) |
||||
|
{ |
||||
|
var tst1 = Math.Abs(matrixH[l - 1, l - 1].Real) + Math.Abs(matrixH[l - 1, l - 1].Imaginary) + Math.Abs(matrixH[l, l].Real) + Math.Abs(matrixH[l, l].Imaginary); |
||||
|
if (Math.Abs(matrixH[l, l - 1].Real) < eps * tst1) |
||||
|
{ |
||||
|
break; |
||||
|
} |
||||
|
|
||||
|
l--; |
||||
|
} |
||||
|
|
||||
|
// Check for convergence
|
||||
|
// One root found
|
||||
|
if (l == n) |
||||
|
{ |
||||
|
matrixH[n, n] += exshift; |
||||
|
vectorV[n] = matrixH[n, n]; |
||||
|
n--; |
||||
|
iter = 0; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
// Form shift
|
||||
|
if (iter != 10 && iter != 20) |
||||
|
{ |
||||
|
s = matrixH[n, n]; |
||||
|
x = matrixH[n - 1, n] * matrixH[n, n - 1].Real; |
||||
|
|
||||
|
if (x.Real != 0.0 || x.Imaginary != 0.0) |
||||
|
{ |
||||
|
y = (matrixH[n - 1, n - 1] - s) / 2.0; |
||||
|
z = ((y * y) + x).SquareRoot(); |
||||
|
if ((y.Real * z.Real) + (y.Imaginary * z.Imaginary) < 0.0) |
||||
|
{ |
||||
|
z *= -1.0; |
||||
|
} |
||||
|
|
||||
|
x /= y + z; |
||||
|
s = s - x; |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
// Form exceptional shift
|
||||
|
s = Math.Abs(matrixH[n, n - 1].Real) + Math.Abs(matrixH[n - 1, n - 2].Real); |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i <= n; i++) |
||||
|
{ |
||||
|
matrixH[i, i] -= s; |
||||
|
} |
||||
|
|
||||
|
exshift += s; |
||||
|
iter++; |
||||
|
|
||||
|
// Reduce to triangle (rows)
|
||||
|
for (var i = l + 1; i <= n; i++) |
||||
|
{ |
||||
|
s = matrixH[i, i - 1].Real; |
||||
|
norm = SpecialFunctions.Hypotenuse(matrixH[i - 1, i - 1].Magnitude, s.Real); |
||||
|
x = matrixH[i - 1, i - 1] / norm; |
||||
|
vectorV[i - 1] = x; |
||||
|
matrixH[i - 1, i - 1] = norm; |
||||
|
matrixH[i, i - 1] = new Complex(0.0, s.Real / norm); |
||||
|
|
||||
|
for (var j = i; j < order; j++) |
||||
|
{ |
||||
|
y = matrixH[i - 1, j]; |
||||
|
z = matrixH[i, j]; |
||||
|
matrixH[i - 1, j] = (x.Conjugate() * y) + (matrixH[i, i - 1].Imaginary * z); |
||||
|
matrixH[i, j] = (x * z) - (matrixH[i, i - 1].Imaginary * y); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
s = matrixH[n, n]; |
||||
|
if (s.Imaginary != 0.0) |
||||
|
{ |
||||
|
s /= matrixH[n, n].Magnitude; |
||||
|
matrixH[n, n] = matrixH[n, n].Magnitude; |
||||
|
|
||||
|
for (var j = n + 1; j < order; j++) |
||||
|
{ |
||||
|
matrixH[n, j] *= s.Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Inverse operation (columns).
|
||||
|
for (var j = l + 1; j <= n; j++) |
||||
|
{ |
||||
|
x = vectorV[j - 1]; |
||||
|
for (var i = 0; i <= j; i++) |
||||
|
{ |
||||
|
z = matrixH[i, j]; |
||||
|
if (i != j) |
||||
|
{ |
||||
|
y = matrixH[i, j - 1]; |
||||
|
matrixH[i, j - 1] = (x * y) + (matrixH[j, j - 1].Imaginary * z); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
y = matrixH[i, j - 1].Real; |
||||
|
matrixH[i, j - 1] = new Complex((x.Real * y.Real) - (x.Imaginary * y.Imaginary) + (matrixH[j, j - 1].Imaginary * z.Real), matrixH[i, j - 1].Imaginary); |
||||
|
} |
||||
|
|
||||
|
matrixH[i, j] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y); |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
y = dataEv[((j - 1) * order) + i]; |
||||
|
z = dataEv[(j * order) + i]; |
||||
|
dataEv[((j - 1) * order) + i] = (x * y) + (matrixH[j, j - 1].Imaginary * z); |
||||
|
dataEv[(j * order) + i] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (s.Imaginary != 0.0) |
||||
|
{ |
||||
|
for (var i = 0; i <= n; i++) |
||||
|
{ |
||||
|
matrixH[i, n] *= s; |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
dataEv[(n * order) + i] *= s; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// All roots found.
|
||||
|
// Backsubstitute to find vectors of upper triangular form
|
||||
|
norm = 0.0; |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
for (var j = i; j < order; j++) |
||||
|
{ |
||||
|
norm = Math.Max(norm, Math.Abs(matrixH[i, j].Real) + Math.Abs(matrixH[i, j].Imaginary)); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (order == 1) |
||||
|
{ |
||||
|
return; |
||||
|
} |
||||
|
|
||||
|
if (norm == 0.0) |
||||
|
{ |
||||
|
return; |
||||
|
} |
||||
|
|
||||
|
for (n = order - 1; n > 0; n--) |
||||
|
{ |
||||
|
x = vectorV[n]; |
||||
|
matrixH[n, n] = 1.0; |
||||
|
|
||||
|
for (var i = n - 1; i >= 0; i--) |
||||
|
{ |
||||
|
z = 0.0; |
||||
|
for (var j = i + 1; j <= n; j++) |
||||
|
{ |
||||
|
z += matrixH[i, j] * matrixH[j, n]; |
||||
|
} |
||||
|
|
||||
|
y = x - vectorV[i]; |
||||
|
if (y.Real == 0.0 && y.Imaginary == 0.0) |
||||
|
{ |
||||
|
y = eps * norm; |
||||
|
} |
||||
|
|
||||
|
matrixH[i, n] = z / y; |
||||
|
|
||||
|
// Overflow control
|
||||
|
var tr = Math.Abs(matrixH[i, n].Real) + Math.Abs(matrixH[i, n].Imaginary); |
||||
|
if ((eps * tr) * tr > 1) |
||||
|
{ |
||||
|
for (var j = i; j <= n; j++) |
||||
|
{ |
||||
|
matrixH[j, n] = matrixH[j, n] / tr; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Back transformation to get eigenvectors of original matrix
|
||||
|
for (var j = order - 1; j > 0; j--) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
z = Complex.Zero; |
||||
|
for (var k = 0; k <= j; k++) |
||||
|
{ |
||||
|
z += dataEv[(k * order) + i] * matrixH[k, j]; |
||||
|
} |
||||
|
|
||||
|
dataEv[(j * order) + i] = z; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
||||
|
public override void Solve(Matrix<Complex> input, Matrix<Complex> 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 (VectorEv.Count != input.RowCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
||||
|
} |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
if (VectorEv.Count != result.RowCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
||||
|
} |
||||
|
|
||||
|
if (IsSymmetric) |
||||
|
{ |
||||
|
var order = VectorEv.Count; |
||||
|
var tmp = new Complex[order]; |
||||
|
|
||||
|
for (var k = 0; k < order; k++) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
Complex value = 0.0; |
||||
|
if (j < order) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(j * order) + i].Conjugate() * input.At(i, k); |
||||
|
} |
||||
|
|
||||
|
value /= VectorEv[j].Real; |
||||
|
} |
||||
|
|
||||
|
tmp[j] = value; |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
Complex value = 0.0; |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(i * order) + j] * tmp[i]; |
||||
|
} |
||||
|
|
||||
|
result[j, k] = value; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSymmetric); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A EVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
||||
|
public override void Solve(Vector<Complex> input, Vector<Complex> result) |
||||
|
{ |
||||
|
if (input == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("input"); |
||||
|
} |
||||
|
|
||||
|
if (result == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("result"); |
||||
|
} |
||||
|
|
||||
|
// Ax=b where A is an m x m matrix
|
||||
|
// Check that b is a column vector with m entries
|
||||
|
if (VectorEv.Count != input.Count) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
||||
|
} |
||||
|
|
||||
|
// Check that x is a column vector with n entries
|
||||
|
if (VectorEv.Count != result.Count) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
||||
|
} |
||||
|
|
||||
|
if (IsSymmetric) |
||||
|
{ |
||||
|
// Symmetric case -> x = V * inv(λ) * VH * b;
|
||||
|
var order = VectorEv.Count; |
||||
|
var tmp = new Complex[order]; |
||||
|
Complex value; |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
value = 0; |
||||
|
if (j < order) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(j * order) + i].Conjugate() * input[i]; |
||||
|
} |
||||
|
|
||||
|
value /= VectorEv[j].Real; |
||||
|
} |
||||
|
|
||||
|
tmp[j] = value; |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
value = 0; |
||||
|
for (int i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(i * order) + j] * tmp[i]; |
||||
|
} |
||||
|
|
||||
|
result[j] = value; |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSymmetric); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Multiply two values T*T
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1">Left operand value</param>
|
||||
|
/// <param name="val2">Right operand value</param>
|
||||
|
/// <returns>Result of multiplication</returns>
|
||||
|
protected sealed override Complex MultiplyT(Complex val1, Complex val2) |
||||
|
{ |
||||
|
return val1 * val2; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,239 @@ |
|||||
|
// <copyright file="DenseGramSchmidt.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using System.Numerics; |
||||
|
using Generic; |
||||
|
using Generic.Factorization; |
||||
|
using Properties; |
||||
|
using Threading; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
||||
|
/// <para>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.</para>
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
||||
|
/// </remarks>
|
||||
|
public class DenseGramSchmidt : GramSchmidt<Complex> |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="DenseGramSchmidt"/> class. This object creates an unitary matrix
|
||||
|
/// using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
|
||||
|
public DenseGramSchmidt(DenseMatrix 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); |
||||
|
Factorize(((DenseMatrix)MatrixQ).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, ((DenseMatrix)MatrixR).Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Factorize matrix using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="q">Initial matrix. On exit is replaced by <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="rowsQ">Number of rows in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="columnsQ">Number of columns in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="r">On exit is filled by <see cref="Matrix{T}"/> R.</param>
|
||||
|
private static void Factorize(Complex[] q, int rowsQ, int columnsQ, Complex[] r) |
||||
|
{ |
||||
|
for (var k = 0; k < columnsQ; k++) |
||||
|
{ |
||||
|
var norm = 0.0; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
norm += q[(k * rowsQ) + i].Magnitude * q[(k * rowsQ) + i].Magnitude; |
||||
|
} |
||||
|
|
||||
|
norm = Math.Sqrt(norm); |
||||
|
if (norm == 0.0) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient); |
||||
|
} |
||||
|
|
||||
|
r[(k * columnsQ) + k] = norm; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
q[(k * rowsQ) + i] /= norm; |
||||
|
} |
||||
|
|
||||
|
for (var j = k + 1; j < columnsQ; j++) |
||||
|
{ |
||||
|
int k1 = k; |
||||
|
int j1 = j; |
||||
|
var dot = CommonParallel.Aggregate(0, rowsQ, index => q[(k1 * rowsQ) + index].Conjugate() * q[(j1 * rowsQ) + index]); |
||||
|
r[(j * columnsQ) + k] = dot; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
var value = q[(j * rowsQ) + i] - (q[(k * rowsQ) + i] * dot); |
||||
|
q[(j * rowsQ) + i] = value; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
||||
|
public override void Solve(Matrix<Complex> input, Matrix<Complex> 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 dinput = input as DenseMatrix; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseMatrix; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
||||
|
public override void Solve(Vector<Complex> input, Vector<Complex> 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 dinput = input as DenseVector; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseVector; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
#region Simple arithmetic of type T
|
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Multiply two values T*T
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1">Left operand value</param>
|
||||
|
/// <param name="val2">Right operand value</param>
|
||||
|
/// <returns>Result of multiplication</returns>
|
||||
|
protected sealed override Complex MultiplyT(Complex val1, Complex val2) |
||||
|
{ |
||||
|
return val1 * val2; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Returns the absolute value of a specified number.
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1"> A number whose absolute is to be found</param>
|
||||
|
/// <returns>Absolute value </returns>
|
||||
|
protected sealed override double AbsoluteT(Complex val1) |
||||
|
{ |
||||
|
return val1.Magnitude; |
||||
|
} |
||||
|
|
||||
|
#endregion
|
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,968 @@ |
|||||
|
// <copyright file="DenseEvd.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using System.Numerics; |
||||
|
using Generic; |
||||
|
using Generic.Factorization; |
||||
|
using Numerics; |
||||
|
using Properties; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Eigenvalues and eigenvectors of a complex matrix.
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// If A is hermitan, then A = V*D*V' where the eigenvalue matrix D is
|
||||
|
/// diagonal and the eigenvector matrix V is hermitan.
|
||||
|
/// I.e. A = V*D*V' and V*VH=I.
|
||||
|
/// If A is not symmetric, then the eigenvalue matrix D is block diagonal
|
||||
|
/// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
||||
|
/// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
|
||||
|
/// columns of V represent the eigenvectors in the sense that A*V = V*D,
|
||||
|
/// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly
|
||||
|
/// conditioned, or even singular, so the validity of the equation
|
||||
|
/// A = V*D*Inverse(V) depends upon V.cond().
|
||||
|
/// </remarks>
|
||||
|
public class DenseEvd : Evd<Complex32> |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="DenseEvd"/> class. This object will compute the
|
||||
|
/// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
|
||||
|
/// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
|
||||
|
public DenseEvd(DenseMatrix matrix) |
||||
|
{ |
||||
|
if (matrix == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("matrix"); |
||||
|
} |
||||
|
|
||||
|
if (matrix.RowCount != matrix.ColumnCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
||||
|
} |
||||
|
|
||||
|
var order = matrix.RowCount; |
||||
|
|
||||
|
// Initialize matricies for eigenvalues and eigenvectors
|
||||
|
MatrixEv = DenseMatrix.Identity(order); |
||||
|
MatrixD = matrix.CreateMatrix(order, order); |
||||
|
VectorEv = new LinearAlgebra.Complex.DenseVector(order); |
||||
|
|
||||
|
IsSymmetric = true; |
||||
|
|
||||
|
for (var i = 0; i < order & IsSymmetric; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < order & IsSymmetric; j++) |
||||
|
{ |
||||
|
IsSymmetric &= matrix[i, j] == matrix[j, i].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (IsSymmetric) |
||||
|
{ |
||||
|
var matrixCopy = matrix.ToArray(); |
||||
|
var tau = new Complex32[order]; |
||||
|
var d = new float[order]; |
||||
|
var e = new float[order]; |
||||
|
|
||||
|
SymmetricTridiagonalize(matrixCopy, d, e, tau, order); |
||||
|
SymmetricDiagonalize(((DenseMatrix)MatrixEv).Data, d, e, order); |
||||
|
SymmetricUntridiagonalize(((DenseMatrix)MatrixEv).Data, matrixCopy, tau, order); |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
VectorEv[i] = new Complex(d[i], e[i]); |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
var matrixH = matrix.ToArray(); |
||||
|
NonsymmetricReduceToHessenberg(((DenseMatrix)MatrixEv).Data, matrixH, order); |
||||
|
NonsymmetricReduceHessenberToRealSchur(((LinearAlgebra.Complex.DenseVector)VectorEv).Data, ((DenseMatrix)MatrixEv).Data, matrixH, order); |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < VectorEv.Count; i++) |
||||
|
{ |
||||
|
MatrixD.At(i, i, (Complex32)VectorEv[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Reduces a complex hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrixA">Source matrix to reduce</param>
|
||||
|
/// <param name="d">Output: Arrays for internal storage of real parts of eigenvalues</param>
|
||||
|
/// <param name="e">Output: Arrays for internal storage of imaginary parts of eigenvalues</param>
|
||||
|
/// <param name="tau">Output: Arrays that contains further information about the transformations.</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures HTRIDI by
|
||||
|
/// Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for
|
||||
|
/// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void SymmetricTridiagonalize(Complex32[,] matrixA, float[] d, float[] e, Complex32[] tau, int order) |
||||
|
{ |
||||
|
float hh; |
||||
|
tau[order - 1] = Complex32.One; |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
d[i] = matrixA[i, i].Real; |
||||
|
} |
||||
|
|
||||
|
// Householder reduction to tridiagonal form.
|
||||
|
for (var i = order - 1; i > 0; i--) |
||||
|
{ |
||||
|
// Scale to avoid under/overflow.
|
||||
|
var scale = 0.0f; |
||||
|
var h = 0.0f; |
||||
|
|
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
scale = scale + Math.Abs(matrixA[i, k].Real) + Math.Abs(matrixA[i, k].Imaginary); |
||||
|
} |
||||
|
|
||||
|
if (scale == 0.0f) |
||||
|
{ |
||||
|
tau[i - 1] = Complex32.One; |
||||
|
e[i] = 0.0f; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
matrixA[i, k] /= scale; |
||||
|
h += matrixA[i, k].MagnitudeSquared; |
||||
|
} |
||||
|
|
||||
|
Complex32 g = (float)Math.Sqrt(h); |
||||
|
e[i] = scale * g.Real; |
||||
|
|
||||
|
Complex32 temp; |
||||
|
var f = matrixA[i, i - 1]; |
||||
|
if (f.Magnitude != 0) |
||||
|
{ |
||||
|
temp = -(matrixA[i, i - 1].Conjugate() * tau[i].Conjugate()) / f.Magnitude; |
||||
|
h += f.Magnitude * g.Real; |
||||
|
g = 1.0f + (g / f.Magnitude); |
||||
|
matrixA[i, i - 1] *= g; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
temp = -tau[i].Conjugate(); |
||||
|
matrixA[i, i - 1] = g; |
||||
|
} |
||||
|
|
||||
|
if ((f.Magnitude == 0) || (i != 1)) |
||||
|
{ |
||||
|
f = Complex32.Zero; |
||||
|
for (var j = 0; j < i; j++) |
||||
|
{ |
||||
|
var tmp = Complex32.Zero; |
||||
|
|
||||
|
// Form element of A*U.
|
||||
|
for (var k = 0; k <= j; k++) |
||||
|
{ |
||||
|
tmp += matrixA[j, k] * matrixA[i, k].Conjugate(); |
||||
|
} |
||||
|
|
||||
|
for (var k = j + 1; k <= i - 1; k++) |
||||
|
{ |
||||
|
tmp += matrixA[k, j].Conjugate() * matrixA[i, k].Conjugate(); |
||||
|
} |
||||
|
|
||||
|
// Form element of P
|
||||
|
tau[j] = tmp / h; |
||||
|
f += (tmp / h) * matrixA[i, j]; |
||||
|
} |
||||
|
|
||||
|
hh = f.Real / (h + h); |
||||
|
|
||||
|
// Form the reduced A.
|
||||
|
for (var j = 0; j < i; j++) |
||||
|
{ |
||||
|
f = matrixA[i, j].Conjugate(); |
||||
|
g = tau[j] - (hh * f); |
||||
|
tau[j] = g.Conjugate(); |
||||
|
|
||||
|
for (var k = 0; k <= j; k++) |
||||
|
{ |
||||
|
matrixA[j, k] -= (f * tau[k]) + (g * matrixA[i, k]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
matrixA[i, k] *= scale; |
||||
|
} |
||||
|
|
||||
|
tau[i - 1] = temp.Conjugate(); |
||||
|
} |
||||
|
|
||||
|
hh = d[i]; |
||||
|
d[i] = matrixA[i, i].Real; |
||||
|
matrixA[i, i] = new Complex32(hh, scale * (float)Math.Sqrt(h)); |
||||
|
} |
||||
|
|
||||
|
hh = d[0]; |
||||
|
d[0] = matrixA[0, 0].Real; |
||||
|
matrixA[0, 0] = hh; |
||||
|
e[0] = 0.0f; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Symmetric tridiagonal QL algorithm.
|
||||
|
/// </summary>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="d">Arrays for internal storage of real parts of eigenvalues</param>
|
||||
|
/// <param name="e">Arrays for internal storage of imaginary parts of eigenvalues</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures tql2, by
|
||||
|
/// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
|
||||
|
/// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void SymmetricDiagonalize(Complex32[] dataEv, float[] d, float[] e, int order) |
||||
|
{ |
||||
|
const int Maxiter = 1000; |
||||
|
|
||||
|
for (var i = 1; i < order; i++) |
||||
|
{ |
||||
|
e[i - 1] = e[i]; |
||||
|
} |
||||
|
|
||||
|
e[order - 1] = 0.0f; |
||||
|
|
||||
|
var f = 0.0f; |
||||
|
var tst1 = 0.0f; |
||||
|
var eps = Precision.DoubleMachinePrecision; |
||||
|
for (var l = 0; l < order; l++) |
||||
|
{ |
||||
|
// Find small subdiagonal element
|
||||
|
tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l])); |
||||
|
var m = l; |
||||
|
while (m < order) |
||||
|
{ |
||||
|
if (Math.Abs(e[m]) <= eps * tst1) |
||||
|
{ |
||||
|
break; |
||||
|
} |
||||
|
|
||||
|
m++; |
||||
|
} |
||||
|
|
||||
|
// If m == l, d[l] is an eigenvalue,
|
||||
|
// otherwise, iterate.
|
||||
|
if (m > l) |
||||
|
{ |
||||
|
var iter = 0; |
||||
|
do |
||||
|
{ |
||||
|
iter = iter + 1; // (Could check iteration count here.)
|
||||
|
|
||||
|
// Compute implicit shift
|
||||
|
var g = d[l]; |
||||
|
var p = (d[l + 1] - g) / (2.0f * e[l]); |
||||
|
var r = SpecialFunctions.Hypotenuse(p, 1.0f); |
||||
|
if (p < 0) |
||||
|
{ |
||||
|
r = -r; |
||||
|
} |
||||
|
|
||||
|
d[l] = e[l] / (p + r); |
||||
|
d[l + 1] = e[l] * (p + r); |
||||
|
|
||||
|
var dl1 = d[l + 1]; |
||||
|
var h = g - d[l]; |
||||
|
for (var i = l + 2; i < order; i++) |
||||
|
{ |
||||
|
d[i] -= h; |
||||
|
} |
||||
|
|
||||
|
f = f + h; |
||||
|
|
||||
|
// Implicit QL transformation.
|
||||
|
p = d[m]; |
||||
|
var c = 1.0f; |
||||
|
var c2 = c; |
||||
|
var c3 = c; |
||||
|
var el1 = e[l + 1]; |
||||
|
var s = 0.0f; |
||||
|
var s2 = 0.0f; |
||||
|
for (var i = m - 1; i >= l; i--) |
||||
|
{ |
||||
|
c3 = c2; |
||||
|
c2 = c; |
||||
|
s2 = s; |
||||
|
g = c * e[i]; |
||||
|
h = c * p; |
||||
|
r = SpecialFunctions.Hypotenuse(p, e[i]); |
||||
|
e[i + 1] = s * r; |
||||
|
s = e[i] / r; |
||||
|
c = p / r; |
||||
|
p = (c * d[i]) - (s * g); |
||||
|
d[i + 1] = h + (s * ((c * g) + (s * d[i]))); |
||||
|
|
||||
|
// Accumulate transformation.
|
||||
|
for (var k = 0; k < order; k++) |
||||
|
{ |
||||
|
h = dataEv[((i + 1) * order) + k].Real; |
||||
|
dataEv[((i + 1) * order) + k] = (s * dataEv[(i * order) + k].Real) + (c * h); |
||||
|
dataEv[(i * order) + k] = (c * dataEv[(i * order) + k].Real) - (s * h); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
p = (-s) * s2 * c3 * el1 * e[l] / dl1; |
||||
|
e[l] = s * p; |
||||
|
d[l] = c * p; |
||||
|
|
||||
|
// Check for convergence. If too many iterations have been performed,
|
||||
|
// throw exception that Convergence Failed
|
||||
|
if (iter >= Maxiter) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ConvergenceFailed); |
||||
|
} |
||||
|
} |
||||
|
while (Math.Abs(e[l]) > eps * tst1); |
||||
|
} |
||||
|
|
||||
|
d[l] = d[l] + f; |
||||
|
e[l] = 0.0f; |
||||
|
} |
||||
|
|
||||
|
// Sort eigenvalues and corresponding vectors.
|
||||
|
for (var i = 0; i < order - 1; i++) |
||||
|
{ |
||||
|
var k = i; |
||||
|
var p = d[i]; |
||||
|
for (var j = i + 1; j < order; j++) |
||||
|
{ |
||||
|
if (d[j] < p) |
||||
|
{ |
||||
|
k = j; |
||||
|
p = d[j]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (k != i) |
||||
|
{ |
||||
|
d[k] = d[i]; |
||||
|
d[i] = p; |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
p = dataEv[(i * order) + j].Real; |
||||
|
dataEv[(i * order) + j] = dataEv[(k * order) + j]; |
||||
|
dataEv[(k * order) + j] = p; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Determines eigenvectors by undoing the symmetric tridiagonalize transformation
|
||||
|
/// </summary>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="matrixA">Previously tridiagonalized matrix by <see cref="SymmetricTridiagonalize"/>.</param>
|
||||
|
/// <param name="tau">Contains further information about the transformations</param>
|
||||
|
/// <param name="order">Input matrix order</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures HTRIBK, by
|
||||
|
/// by Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for
|
||||
|
/// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void SymmetricUntridiagonalize(Complex32[] dataEv, Complex32[,] matrixA, Complex32[] tau, int order) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
dataEv[(j * order) + i] = dataEv[(j * order) + i].Real * tau[i].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Recover and apply the Householder matrices.
|
||||
|
for (var i = 1; i < order; i++) |
||||
|
{ |
||||
|
var h = matrixA[i, i].Imaginary; |
||||
|
if (h != 0) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
var s = Complex32.Zero; |
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
s += dataEv[(j * order) + k] * matrixA[i, k]; |
||||
|
} |
||||
|
|
||||
|
s = (s / h) / h; |
||||
|
|
||||
|
for (var k = 0; k < i; k++) |
||||
|
{ |
||||
|
dataEv[(j * order) + k] -= s * matrixA[i, k].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Nonsymmetric reduction to Hessenberg form.
|
||||
|
/// </summary>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="matrixH">Array for internal storage of nonsymmetric Hessenberg form.</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedures orthes and ortran,
|
||||
|
/// by Martin and Wilkinson, Handbook for Auto. Comp.,
|
||||
|
/// Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutines in EISPACK.</remarks>
|
||||
|
private static void NonsymmetricReduceToHessenberg(Complex32[] dataEv, Complex32[,] matrixH, int order) |
||||
|
{ |
||||
|
var ort = new Complex32[order]; |
||||
|
|
||||
|
for (var m = 1; m < order - 1; m++) |
||||
|
{ |
||||
|
// Scale column.
|
||||
|
var scale = 0.0f; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
scale += Math.Abs(matrixH[i, m - 1].Real) + Math.Abs(matrixH[i, m - 1].Imaginary); |
||||
|
} |
||||
|
|
||||
|
if (scale != 0.0f) |
||||
|
{ |
||||
|
// Compute Householder transformation.
|
||||
|
var h = 0.0f; |
||||
|
for (var i = order - 1; i >= m; i--) |
||||
|
{ |
||||
|
ort[i] = matrixH[i, m - 1] / scale; |
||||
|
h += ort[i].MagnitudeSquared; |
||||
|
} |
||||
|
|
||||
|
var g = (float)Math.Sqrt(h); |
||||
|
if (ort[m].Magnitude != 0) |
||||
|
{ |
||||
|
h = h + (ort[m].Magnitude * g); |
||||
|
g /= ort[m].Magnitude; |
||||
|
ort[m] = (1.0f + g) * ort[m]; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
ort[m] = g; |
||||
|
matrixH[m, m - 1] = scale; |
||||
|
} |
||||
|
|
||||
|
// Apply Householder similarity transformation
|
||||
|
// H = (I-u*u'/h)*H*(I-u*u')/h)
|
||||
|
for (var j = m; j < order; j++) |
||||
|
{ |
||||
|
var f = Complex32.Zero; |
||||
|
for (var i = order - 1; i >= m; i--) |
||||
|
{ |
||||
|
f += ort[i].Conjugate() * matrixH[i, j]; |
||||
|
} |
||||
|
|
||||
|
f = f / h; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
matrixH[i, j] -= f * ort[i]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
var f = Complex32.Zero; |
||||
|
for (var j = order - 1; j >= m; j--) |
||||
|
{ |
||||
|
f += ort[j] * matrixH[i, j]; |
||||
|
} |
||||
|
|
||||
|
f = f / h; |
||||
|
for (var j = m; j < order; j++) |
||||
|
{ |
||||
|
matrixH[i, j] -= f * ort[j].Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
ort[m] = scale * ort[m]; |
||||
|
matrixH[m, m - 1] *= -g; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Accumulate transformations (Algol's ortran).
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
dataEv[(j * order) + i] = i == j ? Complex32.One : Complex32.Zero; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var m = order - 2; m >= 1; m--) |
||||
|
{ |
||||
|
if (matrixH[m, m - 1] != Complex32.Zero && ort[m] != Complex32.Zero) |
||||
|
{ |
||||
|
var norm = (matrixH[m, m - 1].Real * ort[m].Real) + (matrixH[m, m - 1].Imaginary * ort[m].Imaginary); |
||||
|
|
||||
|
for (var i = m + 1; i < order; i++) |
||||
|
{ |
||||
|
ort[i] = matrixH[i, m - 1]; |
||||
|
} |
||||
|
|
||||
|
for (var j = m; j < order; j++) |
||||
|
{ |
||||
|
var g = Complex32.Zero; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
g += ort[i].Conjugate() * dataEv[(j * order) + i]; |
||||
|
} |
||||
|
|
||||
|
// Double division avoids possible underflow
|
||||
|
g /= norm; |
||||
|
for (var i = m; i < order; i++) |
||||
|
{ |
||||
|
dataEv[(j * order) + i] += g * ort[i]; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Create real subdiagonal elements.
|
||||
|
for (var i = 1; i < order; i++) |
||||
|
{ |
||||
|
if (matrixH[i, i - 1].Imaginary != 0.0f) |
||||
|
{ |
||||
|
var y = matrixH[i, i - 1] / matrixH[i, i - 1].Magnitude; |
||||
|
matrixH[i, i - 1] = matrixH[i, i - 1].Magnitude; |
||||
|
for (var j = i; j < order; j++) |
||||
|
{ |
||||
|
matrixH[i, j] *= y.Conjugate(); |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j <= Math.Min(i + 1, order - 1); j++) |
||||
|
{ |
||||
|
matrixH[j, i] *= y; |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
dataEv[(i * order) + j] *= y; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Nonsymmetric reduction from Hessenberg to real Schur form.
|
||||
|
/// </summary>
|
||||
|
/// <param name="vectorV">Data array of the eigenvectors</param>
|
||||
|
/// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
|
||||
|
/// <param name="matrixH">Array for internal storage of nonsymmetric Hessenberg form.</param>
|
||||
|
/// <param name="order">Order of initial matrix</param>
|
||||
|
/// <remarks>This is derived from the Algol procedure hqr2,
|
||||
|
/// by Martin and Wilkinson, Handbook for Auto. Comp.,
|
||||
|
/// Vol.ii-Linear Algebra, and the corresponding
|
||||
|
/// Fortran subroutine in EISPACK.</remarks>
|
||||
|
private static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, Complex32[] dataEv, Complex32[,] matrixH, int order) |
||||
|
{ |
||||
|
// Initialize
|
||||
|
var n = order - 1; |
||||
|
var eps = (float)Precision.SingleMachinePrecision; |
||||
|
|
||||
|
float norm; |
||||
|
Complex32 s, x, y, z, exshift = Complex32.Zero; |
||||
|
|
||||
|
// Outer loop over eigenvalue index
|
||||
|
var iter = 0; |
||||
|
while (n >= 0) |
||||
|
{ |
||||
|
// Look for single small sub-diagonal element
|
||||
|
var l = n; |
||||
|
while (l > 0) |
||||
|
{ |
||||
|
var tst1 = Math.Abs(matrixH[l - 1, l - 1].Real) + Math.Abs(matrixH[l - 1, l - 1].Imaginary) + Math.Abs(matrixH[l, l].Real) + Math.Abs(matrixH[l, l].Imaginary); |
||||
|
if (Math.Abs(matrixH[l, l - 1].Real) < eps * tst1) |
||||
|
{ |
||||
|
break; |
||||
|
} |
||||
|
|
||||
|
l--; |
||||
|
} |
||||
|
|
||||
|
// Check for convergence
|
||||
|
// One root found
|
||||
|
if (l == n) |
||||
|
{ |
||||
|
matrixH[n, n] += exshift; |
||||
|
vectorV[n] = matrixH[n, n].ToComplex(); |
||||
|
n--; |
||||
|
iter = 0; |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
// Form shift
|
||||
|
if (iter != 10 && iter != 20) |
||||
|
{ |
||||
|
s = matrixH[n, n]; |
||||
|
x = matrixH[n - 1, n] * matrixH[n, n - 1].Real; |
||||
|
|
||||
|
if (x.Real != 0.0f || x.Imaginary != 0.0f) |
||||
|
{ |
||||
|
y = (matrixH[n - 1, n - 1] - s) / 2.0f; |
||||
|
z = ((y * y) + x).SquareRoot(); |
||||
|
if ((y.Real * z.Real) + (y.Imaginary * z.Imaginary) < 0.0f) |
||||
|
{ |
||||
|
z *= -1.0f; |
||||
|
} |
||||
|
|
||||
|
x /= y + z; |
||||
|
s = s - x; |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
// Form exceptional shift
|
||||
|
s = Math.Abs(matrixH[n, n - 1].Real) + Math.Abs(matrixH[n - 1, n - 2].Real); |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i <= n; i++) |
||||
|
{ |
||||
|
matrixH[i, i] -= s; |
||||
|
} |
||||
|
|
||||
|
exshift += s; |
||||
|
iter++; |
||||
|
|
||||
|
// Reduce to triangle (rows)
|
||||
|
for (var i = l + 1; i <= n; i++) |
||||
|
{ |
||||
|
s = matrixH[i, i - 1].Real; |
||||
|
norm = SpecialFunctions.Hypotenuse(matrixH[i - 1, i - 1].Magnitude, s.Real); |
||||
|
x = matrixH[i - 1, i - 1] / norm; |
||||
|
vectorV[i - 1] = x.ToComplex(); |
||||
|
matrixH[i - 1, i - 1] = norm; |
||||
|
matrixH[i, i - 1] = new Complex32(0.0f, s.Real / norm); |
||||
|
|
||||
|
for (var j = i; j < order; j++) |
||||
|
{ |
||||
|
y = matrixH[i - 1, j]; |
||||
|
z = matrixH[i, j]; |
||||
|
matrixH[i - 1, j] = (x.Conjugate() * y) + (matrixH[i, i - 1].Imaginary * z); |
||||
|
matrixH[i, j] = (x * z) - (matrixH[i, i - 1].Imaginary * y); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
s = matrixH[n, n]; |
||||
|
if (s.Imaginary != 0.0f) |
||||
|
{ |
||||
|
s /= matrixH[n, n].Magnitude; |
||||
|
matrixH[n, n] = matrixH[n, n].Magnitude; |
||||
|
|
||||
|
for (var j = n + 1; j < order; j++) |
||||
|
{ |
||||
|
matrixH[n, j] *= s.Conjugate(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Inverse operation (columns).
|
||||
|
for (var j = l + 1; j <= n; j++) |
||||
|
{ |
||||
|
x = (Complex32)vectorV[j - 1]; |
||||
|
for (var i = 0; i <= j; i++) |
||||
|
{ |
||||
|
z = matrixH[i, j]; |
||||
|
if (i != j) |
||||
|
{ |
||||
|
y = matrixH[i, j - 1]; |
||||
|
matrixH[i, j - 1] = (x * y) + (matrixH[j, j - 1].Imaginary * z); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
y = matrixH[i, j - 1].Real; |
||||
|
matrixH[i, j - 1] = new Complex32((x.Real * y.Real) - (x.Imaginary * y.Imaginary) + (matrixH[j, j - 1].Imaginary * z.Real), matrixH[i, j - 1].Imaginary); |
||||
|
} |
||||
|
|
||||
|
matrixH[i, j] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y); |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
y = dataEv[((j - 1) * order) + i]; |
||||
|
z = dataEv[(j * order) + i]; |
||||
|
dataEv[((j - 1) * order) + i] = (x * y) + (matrixH[j, j - 1].Imaginary * z); |
||||
|
dataEv[(j * order) + i] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (s.Imaginary != 0.0f) |
||||
|
{ |
||||
|
for (var i = 0; i <= n; i++) |
||||
|
{ |
||||
|
matrixH[i, n] *= s; |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
dataEv[(n * order) + i] *= s; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// All roots found.
|
||||
|
// Backsubstitute to find vectors of upper triangular form
|
||||
|
norm = 0.0f; |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
for (var j = i; j < order; j++) |
||||
|
{ |
||||
|
norm = Math.Max(norm, Math.Abs(matrixH[i, j].Real) + Math.Abs(matrixH[i, j].Imaginary)); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
if (order == 1) |
||||
|
{ |
||||
|
return; |
||||
|
} |
||||
|
|
||||
|
if (norm == 0.0f) |
||||
|
{ |
||||
|
return; |
||||
|
} |
||||
|
|
||||
|
for (n = order - 1; n > 0; n--) |
||||
|
{ |
||||
|
x = (Complex32)vectorV[n]; |
||||
|
matrixH[n, n] = 1.0f; |
||||
|
|
||||
|
for (var i = n - 1; i >= 0; i--) |
||||
|
{ |
||||
|
z = 0.0f; |
||||
|
for (var j = i + 1; j <= n; j++) |
||||
|
{ |
||||
|
z += matrixH[i, j] * matrixH[j, n]; |
||||
|
} |
||||
|
|
||||
|
y = x - (Complex32)vectorV[i]; |
||||
|
if (y.Real == 0.0f && y.Imaginary == 0.0f) |
||||
|
{ |
||||
|
y = eps * norm; |
||||
|
} |
||||
|
|
||||
|
matrixH[i, n] = z / y; |
||||
|
|
||||
|
// Overflow control
|
||||
|
var tr = Math.Abs(matrixH[i, n].Real) + Math.Abs(matrixH[i, n].Imaginary); |
||||
|
if ((eps * tr) * tr > 1) |
||||
|
{ |
||||
|
for (var j = i; j <= n; j++) |
||||
|
{ |
||||
|
matrixH[j, n] = matrixH[j, n] / tr; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Back transformation to get eigenvectors of original matrix
|
||||
|
for (var j = order - 1; j > 0; j--) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
z = Complex32.Zero; |
||||
|
for (var k = 0; k <= j; k++) |
||||
|
{ |
||||
|
z += dataEv[(k * order) + i] * matrixH[k, j]; |
||||
|
} |
||||
|
|
||||
|
dataEv[(j * order) + i] = z; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
||||
|
public override void Solve(Matrix<Complex32> input, Matrix<Complex32> 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 (VectorEv.Count != input.RowCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
||||
|
} |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
if (VectorEv.Count != result.RowCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
||||
|
} |
||||
|
|
||||
|
if (IsSymmetric) |
||||
|
{ |
||||
|
var order = VectorEv.Count; |
||||
|
var tmp = new Complex32[order]; |
||||
|
|
||||
|
for (var k = 0; k < order; k++) |
||||
|
{ |
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
Complex32 value = 0.0f; |
||||
|
if (j < order) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(j * order) + i].Conjugate() * input.At(i, k); |
||||
|
} |
||||
|
|
||||
|
value /= (float)VectorEv[j].Real; |
||||
|
} |
||||
|
|
||||
|
tmp[j] = value; |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
Complex32 value = 0.0f; |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(i * order) + j] * tmp[i]; |
||||
|
} |
||||
|
|
||||
|
result[j, k] = value; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSymmetric); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A EVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
||||
|
public override void Solve(Vector<Complex32> input, Vector<Complex32> result) |
||||
|
{ |
||||
|
if (input == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("input"); |
||||
|
} |
||||
|
|
||||
|
if (result == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("result"); |
||||
|
} |
||||
|
|
||||
|
// Ax=b where A is an m x m matrix
|
||||
|
// Check that b is a column vector with m entries
|
||||
|
if (VectorEv.Count != input.Count) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
||||
|
} |
||||
|
|
||||
|
// Check that x is a column vector with n entries
|
||||
|
if (VectorEv.Count != result.Count) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
||||
|
} |
||||
|
|
||||
|
if (IsSymmetric) |
||||
|
{ |
||||
|
// Symmetric case -> x = V * inv(λ) * VH * b;
|
||||
|
var order = VectorEv.Count; |
||||
|
var tmp = new Complex32[order]; |
||||
|
Complex32 value; |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
value = 0; |
||||
|
if (j < order) |
||||
|
{ |
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(j * order) + i].Conjugate() * input[i]; |
||||
|
} |
||||
|
|
||||
|
value /= (float)VectorEv[j].Real; |
||||
|
} |
||||
|
|
||||
|
tmp[j] = value; |
||||
|
} |
||||
|
|
||||
|
for (var j = 0; j < order; j++) |
||||
|
{ |
||||
|
value = 0; |
||||
|
for (int i = 0; i < order; i++) |
||||
|
{ |
||||
|
value += ((DenseMatrix)MatrixEv).Data[(i * order) + j] * tmp[i]; |
||||
|
} |
||||
|
|
||||
|
result[j] = value; |
||||
|
} |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSymmetric); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Multiply two values T*T
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1">Left operand value</param>
|
||||
|
/// <param name="val2">Right operand value</param>
|
||||
|
/// <returns>Result of multiplication</returns>
|
||||
|
protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) |
||||
|
{ |
||||
|
return val1 * val2; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,239 @@ |
|||||
|
// <copyright file="DenseGramSchmidt.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using Generic; |
||||
|
using Generic.Factorization; |
||||
|
using Numerics; |
||||
|
using Properties; |
||||
|
using Threading; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
||||
|
/// <para>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.</para>
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
||||
|
/// </remarks>
|
||||
|
public class DenseGramSchmidt : GramSchmidt<Complex32> |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="DenseGramSchmidt"/> class. This object creates an unitary matrix
|
||||
|
/// using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
|
||||
|
public DenseGramSchmidt(DenseMatrix 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); |
||||
|
Factorize(((DenseMatrix)MatrixQ).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, ((DenseMatrix)MatrixR).Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Factorize matrix using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="q">Initial matrix. On exit is replaced by <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="rowsQ">Number of rows in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="columnsQ">Number of columns in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="r">On exit is filled by <see cref="Matrix{T}"/> R.</param>
|
||||
|
private static void Factorize(Complex32[] q, int rowsQ, int columnsQ, Complex32[] r) |
||||
|
{ |
||||
|
for (var k = 0; k < columnsQ; k++) |
||||
|
{ |
||||
|
var norm = 0.0f; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
norm += q[(k * rowsQ) + i].Magnitude * q[(k * rowsQ) + i].Magnitude; |
||||
|
} |
||||
|
|
||||
|
norm = (float)Math.Sqrt(norm); |
||||
|
if (norm == 0.0f) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient); |
||||
|
} |
||||
|
|
||||
|
r[(k * columnsQ) + k] = norm; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
q[(k * rowsQ) + i] /= norm; |
||||
|
} |
||||
|
|
||||
|
for (var j = k + 1; j < columnsQ; j++) |
||||
|
{ |
||||
|
int k1 = k; |
||||
|
int j1 = j; |
||||
|
var dot = CommonParallel.Aggregate(0, rowsQ, index => q[(k1 * rowsQ) + index].Conjugate() * q[(j1 * rowsQ) + index]); |
||||
|
r[(j * columnsQ) + k] = dot; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
var value = q[(j * rowsQ) + i] - (q[(k * rowsQ) + i] * dot); |
||||
|
q[(j * rowsQ) + i] = value; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
||||
|
public override void Solve(Matrix<Complex32> input, Matrix<Complex32> 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 dinput = input as DenseMatrix; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseMatrix; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
||||
|
public override void Solve(Vector<Complex32> input, Vector<Complex32> 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 dinput = input as DenseVector; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseVector; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
#region Simple arithmetic of type T
|
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Multiply two values T*T
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1">Left operand value</param>
|
||||
|
/// <param name="val2">Right operand value</param>
|
||||
|
/// <returns>Result of multiplication</returns>
|
||||
|
protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) |
||||
|
{ |
||||
|
return val1 * val2; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Returns the absolute value of a specified number.
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1"> A number whose absolute is to be found</param>
|
||||
|
/// <returns>Absolute value </returns>
|
||||
|
protected sealed override double AbsoluteT(Complex32 val1) |
||||
|
{ |
||||
|
return val1.Magnitude; |
||||
|
} |
||||
|
|
||||
|
#endregion
|
||||
|
} |
||||
|
} |
||||
File diff suppressed because it is too large
@ -0,0 +1,237 @@ |
|||||
|
// <copyright file="DenseGramSchmidt.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using Generic; |
||||
|
using Generic.Factorization; |
||||
|
using Properties; |
||||
|
using Threading; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
||||
|
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix.</para>
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
||||
|
/// </remarks>
|
||||
|
public class DenseGramSchmidt : GramSchmidt<double> |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="DenseGramSchmidt"/> class. This object creates an orthogonal matrix
|
||||
|
/// using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
|
||||
|
public DenseGramSchmidt(DenseMatrix 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); |
||||
|
Factorize(((DenseMatrix)MatrixQ).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, ((DenseMatrix)MatrixR).Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Factorize matrix using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="q">Initial matrix. On exit is replaced by <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="rowsQ">Number of rows in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="columnsQ">Number of columns in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="r">On exit is filled by <see cref="Matrix{T}"/> R.</param>
|
||||
|
private static void Factorize(double[] q, int rowsQ, int columnsQ, double[] r) |
||||
|
{ |
||||
|
for (var k = 0; k < columnsQ; k++) |
||||
|
{ |
||||
|
var norm = 0.0; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
norm += q[(k * rowsQ) + i] * q[(k * rowsQ) + i]; |
||||
|
} |
||||
|
|
||||
|
norm = Math.Sqrt(norm); |
||||
|
if (norm == 0.0) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient); |
||||
|
} |
||||
|
|
||||
|
r[(k * columnsQ) + k] = norm; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
q[(k * rowsQ) + i] /= norm; |
||||
|
} |
||||
|
|
||||
|
for (var j = k + 1; j < columnsQ; j++) |
||||
|
{ |
||||
|
int k1 = k; |
||||
|
int j1 = j; |
||||
|
var dot = CommonParallel.Aggregate(0, rowsQ, index => q[(k1 * rowsQ) + index] * q[(j1 * rowsQ) + index]); |
||||
|
r[(j * columnsQ) + k] = dot; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
var value = q[(j * rowsQ) + i] - (q[(k * rowsQ) + i] * dot); |
||||
|
q[(j * rowsQ) + i] = value; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
||||
|
public override void Solve(Matrix<double> input, Matrix<double> 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 dinput = input as DenseMatrix; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseMatrix; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
||||
|
public override void Solve(Vector<double> input, Vector<double> 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 dinput = input as DenseVector; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseVector; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
#region Simple arithmetic of type T
|
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Multiply two values T*T
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1">Left operand value</param>
|
||||
|
/// <param name="val2">Right operand value</param>
|
||||
|
/// <returns>Result of multiplication</returns>
|
||||
|
protected sealed override double MultiplyT(double val1, double val2) |
||||
|
{ |
||||
|
return val1 * val2; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Returns the absolute value of a specified number.
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1"> A number whose absolute is to be found</param>
|
||||
|
/// <returns>Absolute value </returns>
|
||||
|
protected sealed override double AbsoluteT(double val1) |
||||
|
{ |
||||
|
return Math.Abs(val1); |
||||
|
} |
||||
|
#endregion
|
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,138 @@ |
|||||
|
// <copyright file="GramSchmidt.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using System.Numerics; |
||||
|
using Generic; |
||||
|
using Numerics; |
||||
|
using Properties; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
||||
|
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix.</para>
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
||||
|
/// </remarks>
|
||||
|
/// <typeparam name="T">Supported data types are double, single, <see cref="Complex"/>, and <see cref="Complex32"/>.</typeparam>
|
||||
|
public abstract class GramSchmidt<T> : QR<T> |
||||
|
where T : struct, IEquatable<T>, IFormattable |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Internal method which routes the call to perform the QR factorization to the appropriate class.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <returns>A QR factorization object.</returns>
|
||||
|
new internal static GramSchmidt<T> Create(Matrix<T> matrix) |
||||
|
{ |
||||
|
if (typeof(T) == typeof(double)) |
||||
|
{ |
||||
|
var dense = matrix as LinearAlgebra.Double.DenseMatrix; |
||||
|
if (dense != null) |
||||
|
{ |
||||
|
return new LinearAlgebra.Double.Factorization.DenseGramSchmidt(dense) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
return new LinearAlgebra.Double.Factorization.UserGramSchmidt(matrix as Matrix<double>) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
if (typeof(T) == typeof(float)) |
||||
|
{ |
||||
|
var dense = matrix as LinearAlgebra.Single.DenseMatrix; |
||||
|
if (dense != null) |
||||
|
{ |
||||
|
return new LinearAlgebra.Single.Factorization.DenseGramSchmidt(dense) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
return new LinearAlgebra.Single.Factorization.UserGramSchmidt(matrix as Matrix<float>) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
if (typeof(T) == typeof(Complex)) |
||||
|
{ |
||||
|
var dense = matrix as LinearAlgebra.Complex.DenseMatrix; |
||||
|
if (dense != null) |
||||
|
{ |
||||
|
return new LinearAlgebra.Complex.Factorization.DenseGramSchmidt(dense) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
return new LinearAlgebra.Complex.Factorization.UserGramSchmidt(matrix as Matrix<Complex>) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
if (typeof(T) == typeof(Complex32)) |
||||
|
{ |
||||
|
var dense = matrix as LinearAlgebra.Complex32.DenseMatrix; |
||||
|
if (dense != null) |
||||
|
{ |
||||
|
return new LinearAlgebra.Complex32.Factorization.DenseGramSchmidt(dense) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
return new LinearAlgebra.Complex32.Factorization.UserGramSchmidt(matrix as Matrix<Complex32>) as GramSchmidt<T>; |
||||
|
} |
||||
|
|
||||
|
throw new NotImplementedException(); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets a value indicating whether the matrix is full rank or not.
|
||||
|
/// </summary>
|
||||
|
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
||||
|
public sealed override bool IsFullRank |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return true; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
||||
|
/// </summary>
|
||||
|
public override double Determinant |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
if (MatrixQ.RowCount != MatrixQ.ColumnCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
||||
|
} |
||||
|
|
||||
|
var det = OneValueT; |
||||
|
for (var i = 0; i < MatrixR.ColumnCount; i++) |
||||
|
{ |
||||
|
det = MultiplyT(det, MatrixR.At(i, i)); |
||||
|
if (AbsoluteT(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, (typeof(T) == typeof(float) || typeof(T) == typeof(Complex32)) ? 7 : 15)) |
||||
|
{ |
||||
|
return 0; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
return AbsoluteT(det); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
File diff suppressed because it is too large
@ -0,0 +1,238 @@ |
|||||
|
// <copyright file="DenseGramSchmidt.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using Generic; |
||||
|
using Generic.Factorization; |
||||
|
using Properties; |
||||
|
using Threading; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
||||
|
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix.</para>
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
||||
|
/// </remarks>
|
||||
|
public class DenseGramSchmidt : GramSchmidt<float> |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="DenseGramSchmidt"/> class. This object creates an orthogonal matrix
|
||||
|
/// using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
|
||||
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is rank deficient</exception>
|
||||
|
public DenseGramSchmidt(DenseMatrix 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); |
||||
|
Factorize(((DenseMatrix)MatrixQ).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, ((DenseMatrix)MatrixR).Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Factorize matrix using the modified Gram-Schmidt method.
|
||||
|
/// </summary>
|
||||
|
/// <param name="q">Initial matrix. On exit is replaced by <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="rowsQ">Number of rows in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="columnsQ">Number of columns in <see cref="Matrix{T}"/> Q.</param>
|
||||
|
/// <param name="r">On exit is filled by <see cref="Matrix{T}"/> R.</param>
|
||||
|
private static void Factorize(float[] q, int rowsQ, int columnsQ, float[] r) |
||||
|
{ |
||||
|
for (var k = 0; k < columnsQ; k++) |
||||
|
{ |
||||
|
var norm = 0.0f; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
norm += q[(k * rowsQ) + i] * q[(k * rowsQ) + i]; |
||||
|
} |
||||
|
|
||||
|
norm = (float)Math.Sqrt(norm); |
||||
|
if (norm == 0.0) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient); |
||||
|
} |
||||
|
|
||||
|
r[(k * columnsQ) + k] = norm; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
q[(k * rowsQ) + i] /= norm; |
||||
|
} |
||||
|
|
||||
|
for (var j = k + 1; j < columnsQ; j++) |
||||
|
{ |
||||
|
int k1 = k; |
||||
|
int j1 = j; |
||||
|
var dot = CommonParallel.Aggregate(0, rowsQ, index => q[(k1 * rowsQ) + index] * q[(j1 * rowsQ) + index]); |
||||
|
r[(j * columnsQ) + k] = dot; |
||||
|
for (var i = 0; i < rowsQ; i++) |
||||
|
{ |
||||
|
var value = q[(j * rowsQ) + i] - (q[(k * rowsQ) + i] * dot); |
||||
|
q[(j * rowsQ) + i] = value; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
||||
|
public override void Solve(Matrix<float> input, Matrix<float> 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 dinput = input as DenseMatrix; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseMatrix; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A QR factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
||||
|
public override void Solve(Vector<float> input, Vector<float> 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 dinput = input as DenseVector; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseVector; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do GramSchmidt factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
#region Simple arithmetic of type T
|
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Multiply two values T*T
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1">Left operand value</param>
|
||||
|
/// <param name="val2">Right operand value</param>
|
||||
|
/// <returns>Result of multiplication</returns>
|
||||
|
protected sealed override float MultiplyT(float val1, float val2) |
||||
|
{ |
||||
|
return val1 * val2; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Returns the absolute value of a specified number.
|
||||
|
/// </summary>
|
||||
|
/// <param name="val1"> A number whose absolute is to be found</param>
|
||||
|
/// <returns>Absolute value </returns>
|
||||
|
protected sealed override double AbsoluteT(float val1) |
||||
|
{ |
||||
|
return Math.Abs(val1); |
||||
|
} |
||||
|
|
||||
|
#endregion
|
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,219 @@ |
|||||
|
// <copyright file="PoissonTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.DistributionTests.Discrete |
||||
|
{ |
||||
|
using System; |
||||
|
using System.Linq; |
||||
|
using MbUnit.Framework; |
||||
|
using Distributions; |
||||
|
|
||||
|
[TestFixture] |
||||
|
public class PoissonTests |
||||
|
{ |
||||
|
[SetUp] |
||||
|
public void SetUp() |
||||
|
{ |
||||
|
Control.CheckDistributionParameters = true; |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5)] |
||||
|
[Row(5.4)] |
||||
|
[Row(10.8)] |
||||
|
public void CanCreatePoisson(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreEqual(lambda, d.Lambda); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedException(typeof(ArgumentOutOfRangeException))] |
||||
|
[Row(Double.NaN)] |
||||
|
[Row(-1.5)] |
||||
|
[Row(0.0)] |
||||
|
public void PoissonCreateFailsWithBadParameters(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
public void ValidateToString() |
||||
|
{ |
||||
|
var d = new Poisson(0.3); |
||||
|
Assert.AreEqual(String.Format("Poisson(λ = {0})", 0.3), d.ToString()); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5)] |
||||
|
[Row(5.4)] |
||||
|
[Row(10.8)] |
||||
|
public void CanSetProbabilityOfOne(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(0.3) |
||||
|
{ |
||||
|
Lambda = lambda |
||||
|
}; |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedException(typeof(ArgumentOutOfRangeException))] |
||||
|
[Row(Double.NaN)] |
||||
|
[Row(-1.5)] |
||||
|
[Row(0.0)] |
||||
|
public void SetProbabilityOfOneFails(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(0.3) |
||||
|
{ |
||||
|
Lambda = lambda |
||||
|
}; |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5)] |
||||
|
[Row(5.4)] |
||||
|
[Row(10.8)] |
||||
|
public void ValidateEntropy(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreEqual((0.5 * Math.Log(2 * Math.PI * Math.E * lambda)) - (1.0 / (12.0 * lambda)) - 1.0 / (24.0 * lambda * lambda) - (19.0 / (360.0 * lambda * lambda * lambda)), d.Entropy); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5)] |
||||
|
[Row(5.4)] |
||||
|
[Row(10.8)] |
||||
|
public void ValidateSkewness(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreEqual(1.0 / Math.Sqrt(lambda), d.Skewness); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5)] |
||||
|
[Row(5.4)] |
||||
|
[Row(10.8)] |
||||
|
public void ValidateMode(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreEqual((int)Math.Floor(lambda), d.Mode); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5)] |
||||
|
[Row(5.4)] |
||||
|
[Row(10.8)] |
||||
|
public void ValidateMedian(double lambda) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreEqual((int)Math.Floor(lambda + (1.0 / 3.0) - (0.02 / lambda)), d.Median); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
public void ValidateMinimum() |
||||
|
{ |
||||
|
var d = new Poisson(0.3); |
||||
|
Assert.AreEqual(0, d.Minimum); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
public void ValidateMaximum() |
||||
|
{ |
||||
|
var d = new Poisson(0.3); |
||||
|
Assert.AreEqual(int.MaxValue, d.Maximum); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5, 1, 0.334695240222645000000000000000)] |
||||
|
[Row(1.5, 10, 0.000003545747740570180000000000)] |
||||
|
[Row(1.5, 20, 0.000000000000000304971208961018)] |
||||
|
[Row(5.4, 1, 0.024389537090108400000000000000)] |
||||
|
[Row(5.4, 10, 0.026241240591792300000000000000)] |
||||
|
[Row(5.4, 20, 0.000000825202200316548000000000)] |
||||
|
[Row(10.8, 1, 0.000220314636840657000000000000)] |
||||
|
[Row(10.8, 10, 0.121365183659420000000000000000)] |
||||
|
[Row(10.8, 20, 0.003908139778574110000000000000)] |
||||
|
public void ValidateProbability(double lambda, int x, double result) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreApproximatelyEqual(d.Probability(x), result, 1e-12); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5, 1, 0.334695240222645000000000000000)] |
||||
|
[Row(1.5, 10, 0.000003545747740570180000000000)] |
||||
|
[Row(1.5, 20, 0.000000000000000304971208961018)] |
||||
|
[Row(5.4, 1, 0.024389537090108400000000000000)] |
||||
|
[Row(5.4, 10, 0.026241240591792300000000000000)] |
||||
|
[Row(5.4, 20, 0.000000825202200316548000000000)] |
||||
|
[Row(10.8, 1, 0.000220314636840657000000000000)] |
||||
|
[Row(10.8, 10, 0.121365183659420000000000000000)] |
||||
|
[Row(10.8, 20, 0.003908139778574110000000000000)] |
||||
|
public void ValidateProbabilityLn(double lambda, int x, double result) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreApproximatelyEqual(d.ProbabilityLn(x), Math.Log(result), 1e-12); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
public void CanSample() |
||||
|
{ |
||||
|
var d = new Poisson(0.3); |
||||
|
var s = d.Sample(); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
public void CanSampleSequence() |
||||
|
{ |
||||
|
var d = new Poisson(0.3); |
||||
|
var ied = d.Samples(); |
||||
|
var e = ied.Take(5).ToArray(); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1.5, 1, 0.5578254003710750000000)] |
||||
|
[Row(1.5, 10, 0.9999994482467640000000)] |
||||
|
[Row(1.5, 20, 1.0000000000000000000000)] |
||||
|
[Row(5.4, 1, 0.0289061180327211000000)] |
||||
|
[Row(5.4, 10, 0.9774863006897650000000)] |
||||
|
[Row(5.4, 20, 0.9999997199928290000000)] |
||||
|
[Row(10.8, 1, 0.0002407141402518290000)] |
||||
|
[Row(10.8, 10, 0.4839692359955690000000)] |
||||
|
[Row(10.8, 20, 0.9961800769608090000000)] |
||||
|
[Row(20.1, 1, 0.0000000393516882521484)] |
||||
|
[Row(20.1, 10, 0.0102444128791257000000)] |
||||
|
[Row(20.1, 20, 0.5502097908860160000000)] |
||||
|
public void ValidateCumulativeDistribution(double lambda, int x, double result) |
||||
|
{ |
||||
|
var d = new Poisson(lambda); |
||||
|
Assert.AreApproximatelyEqual(d.CumulativeDistribution(x), result, 1e-12); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,364 @@ |
|||||
|
// <copyright file="EvdTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization |
||||
|
{ |
||||
|
using System.Numerics; |
||||
|
using LinearAlgebra.Complex; |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Complex.Factorization; |
||||
|
|
||||
|
public class EvdTests |
||||
|
{ |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new DenseEvd(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorEvd.EValues().Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A*V = λ*V
|
||||
|
var matrixAv = matrixA * factorEvd.EVectors(); |
||||
|
var matrixLv = factorEvd.EVectors() * factorEvd.D(); |
||||
|
|
||||
|
for (var i = 0; i < matrixAv.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixAv.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixAv[i, j].Real, matrixLv[i, j].Real, 1e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixAv[i, j].Imaginary, matrixLv[i, j].Imaginary, 1e-9); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A = V*λ*VT
|
||||
|
var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().ConjugateTranspose(); |
||||
|
|
||||
|
for (var i = 0; i < matrix.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrix.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrix[i, j].Real, matrixA[i, j].Real, 1e-9); |
||||
|
Assert.AreApproximatelyEqual(matrix[i, j].Imaginary, matrixA[i, j].Imaginary, 1e-9); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankSquare(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Rank, order); |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankOfSquareSingular(int order) |
||||
|
{ |
||||
|
var matrixA = new DenseMatrix(order, order); |
||||
|
matrixA[0, 0] = 1; |
||||
|
matrixA[order - 1, order - 1] = 1; |
||||
|
for (var i = 1; i < order - 1; i++) |
||||
|
{ |
||||
|
matrixA[i, i - 1] = 1; |
||||
|
matrixA[i, i + 1] = 1; |
||||
|
matrixA[i - 1, i] = 1; |
||||
|
matrixA[i + 1, i] = 1; |
||||
|
} |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Determinant, 0); |
||||
|
Assert.AreEqual(factorEvd.Rank, order - 1); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
Assert.AreEqual(1.0, factorEvd.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var resultx = factorEvd.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1e-9); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1e-9); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixX = factorEvd.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1e-9); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrixWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new DenseVector(order); |
||||
|
factorEvd.Solve(vectorb, resultx); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1e-9); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1e-9); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new DenseMatrix(order, order); |
||||
|
factorEvd.Solve(matrixB, matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1e-9); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,356 @@ |
|||||
|
// <copyright file="UserGramSchmidtTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization |
||||
|
{ |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Complex.Factorization; |
||||
|
|
||||
|
public class UserGramSchmidtTests |
||||
|
{ |
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new UserGramSchmidt(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentException] |
||||
|
public void WideMatrixThrowsInvalidMatrixOperationException() |
||||
|
{ |
||||
|
new UserGramSchmidt(new UserDefinedMatrix(3, 4)); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.Q.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.Q.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1,1)] |
||||
|
[Row(2,2)] |
||||
|
[Row(5,5)] |
||||
|
[Row(10,6)] |
||||
|
[Row(50,48)] |
||||
|
[Row(100,98)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
// Make sure the Q has the right dimensions.
|
||||
|
Assert.AreEqual(row, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
// Make sure the R has the right dimensions.
|
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.ColumnCount); |
||||
|
|
||||
|
// Make sure the R factor is upper triangular.
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i > j) |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure the Q*R is the original matrix.
|
||||
|
var matrixQfromR = factorGramSchmidt.Q * factorGramSchmidt.R; |
||||
|
for (var i = 0; i < matrixQfromR.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixQfromR.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixA[i, j].Real, matrixQfromR[i, j].Real, 1.0e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixA[i, j].Imaginary, matrixQfromR[i, j].Imaginary, 1.0e-9); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure the Q is unitary --> (Q*)x(Q) = I
|
||||
|
var matrixQсtQ = factorGramSchmidt.Q.ConjugateTranspose() * factorGramSchmidt.Q; |
||||
|
for (var i = 0; i < matrixQсtQ.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixQсtQ.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Real, 1.0, 1.0e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Imaginary, 0.0, 1.0e-9); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Real, 0.0, 1.0e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Imaginary, 0.0, 1.0e-9); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVector(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var resultx = factorGramSchmidt.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1.0e-9); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1.0e-9); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixX = factorGramSchmidt.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1.0e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1.0e-9); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new UserDefinedVector(order); |
||||
|
factorGramSchmidt.Solve(vectorb,resultx); |
||||
|
|
||||
|
Assert.AreEqual(vectorb.Count, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1.0e-9); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1.0e-9); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new UserDefinedMatrix(order, order); |
||||
|
factorGramSchmidt.Solve(matrixB,matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1.0e-9); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1.0e-9); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,365 @@ |
|||||
|
// <copyright file="EvdTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization |
||||
|
{ |
||||
|
using System.Numerics; |
||||
|
using LinearAlgebra.Complex32; |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Complex.Factorization; |
||||
|
|
||||
|
public class EvdTests |
||||
|
{ |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new DenseEvd(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorEvd.EValues().Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A*V = λ*V
|
||||
|
var matrixAv = matrixA * factorEvd.EVectors(); |
||||
|
var matrixLv = factorEvd.EVectors() * factorEvd.D(); |
||||
|
|
||||
|
for (var i = 0; i < matrixAv.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixAv.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixAv[i, j].Real, matrixLv[i, j].Real, 1e-4f); |
||||
|
Assert.AreApproximatelyEqual(matrixAv[i, j].Imaginary, matrixLv[i, j].Imaginary, 1e-4f); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
[Ignore] |
||||
|
public void CanFactorizeRandomSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A = V*λ*VT
|
||||
|
var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().ConjugateTranspose(); |
||||
|
|
||||
|
for (var i = 0; i < matrix.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrix.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrix[i, j].Real, matrixA[i, j].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(matrix[i, j].Imaginary, matrixA[i, j].Imaginary, 1e-3f); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankSquare(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Rank, order); |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankOfSquareSingular(int order) |
||||
|
{ |
||||
|
var matrixA = new DenseMatrix(order, order); |
||||
|
matrixA[0, 0] = 1; |
||||
|
matrixA[order - 1, order - 1] = 1; |
||||
|
for (var i = 1; i < order - 1; i++) |
||||
|
{ |
||||
|
matrixA[i, i - 1] = 1; |
||||
|
matrixA[i, i + 1] = 1; |
||||
|
matrixA[i - 1, i] = 1; |
||||
|
matrixA[i + 1, i] = 1; |
||||
|
} |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Determinant, 0); |
||||
|
Assert.AreEqual(factorEvd.Rank, order - 1); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
Assert.AreEqual(1.0, factorEvd.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
//[Row(100)]
|
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var resultx = factorEvd.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1e-3f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixX = factorEvd.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1e-2f); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1e-2f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
// [Row(100)]
|
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrixWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new DenseVector(order); |
||||
|
factorEvd.Solve(vectorb, resultx); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1e-3f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new DenseMatrix(order, order); |
||||
|
factorEvd.Solve(matrixB, matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1e-2f); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1e-2f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,356 @@ |
|||||
|
// <copyright file="UserGramSchmidtTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization |
||||
|
{ |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Complex32.Factorization; |
||||
|
|
||||
|
public class UserGramSchmidtTests |
||||
|
{ |
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new UserGramSchmidt(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentException] |
||||
|
public void WideMatrixThrowsInvalidMatrixOperationException() |
||||
|
{ |
||||
|
new UserGramSchmidt(new UserDefinedMatrix(3, 4)); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0f, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0f, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.Q.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.Q.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0f, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0f, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
Assert.AreEqual(1.0f, factorGramSchmidt.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1,1)] |
||||
|
[Row(2,2)] |
||||
|
[Row(5,5)] |
||||
|
[Row(10,6)] |
||||
|
[Row(50,48)] |
||||
|
[Row(100,98)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
// Make sure the Q has the right dimensions.
|
||||
|
Assert.AreEqual(row, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
// Make sure the R has the right dimensions.
|
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.ColumnCount); |
||||
|
|
||||
|
// Make sure the R factor is upper triangular.
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i > j) |
||||
|
{ |
||||
|
Assert.AreEqual(0.0f, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure the Q*R is the original matrix.
|
||||
|
var matrixQfromR = factorGramSchmidt.Q * factorGramSchmidt.R; |
||||
|
for (var i = 0; i < matrixQfromR.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixQfromR.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixA[i, j].Real, matrixQfromR[i, j].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(matrixA[i, j].Imaginary, matrixQfromR[i, j].Imaginary, 1e-3f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure the Q is unitary --> (Q*)x(Q) = I
|
||||
|
var matrixQсtQ = factorGramSchmidt.Q.ConjugateTranspose() * factorGramSchmidt.Q; |
||||
|
for (var i = 0; i < matrixQсtQ.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixQсtQ.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Real, 1.0f, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Imaginary, 0.0f, 1e-3f); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Real, 0.0f, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(matrixQсtQ[i, j].Imaginary, 0.0f, 1e-3f); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVector(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var resultx = factorGramSchmidt.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1e-3f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixX = factorGramSchmidt.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1e-3f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new UserDefinedVector(order); |
||||
|
factorGramSchmidt.Solve(vectorb,resultx); |
||||
|
|
||||
|
Assert.AreEqual(vectorb.Count, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Real, bReconstruct[i].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(vectorb[i].Imaginary, bReconstruct[i].Imaginary, 1e-3f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new UserDefinedMatrix(order, order); |
||||
|
factorGramSchmidt.Solve(matrixB,matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Real, matrixBReconstruct[i, j].Real, 1e-3f); |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j].Imaginary, matrixBReconstruct[i, j].Imaginary, 1e-3f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,358 @@ |
|||||
|
// <copyright file="EvdTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization |
||||
|
{ |
||||
|
using System.Numerics; |
||||
|
using LinearAlgebra.Double; |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Double.Factorization; |
||||
|
|
||||
|
public class EvdTests |
||||
|
{ |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new DenseEvd(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorEvd.EValues().Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A*V = λ*V
|
||||
|
var matrixAv = matrixA * factorEvd.EVectors(); |
||||
|
var matrixLv = factorEvd.EVectors() * factorEvd.D(); |
||||
|
|
||||
|
for (var i = 0; i < matrixAv.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixAv.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixAv[i, j], matrixLv[i, j], 1.0e-10); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A = V*λ*VT
|
||||
|
var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().Transpose(); |
||||
|
|
||||
|
for (var i = 0; i < matrix.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrix.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrix[i, j], matrixA[i, j], 1.0e-10); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankSquare(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Rank, order); |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankOfSquareSingular(int order) |
||||
|
{ |
||||
|
var matrixA = new DenseMatrix(order, order); |
||||
|
matrixA[0, 0] = 1; |
||||
|
matrixA[order - 1, order - 1] = 1; |
||||
|
for (var i = 1; i < order - 1; i++) |
||||
|
{ |
||||
|
matrixA[i, i - 1] = 1; |
||||
|
matrixA[i, i + 1] = 1; |
||||
|
matrixA[i - 1, i] = 1; |
||||
|
matrixA[i + 1, i] = 1; |
||||
|
} |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Determinant, 0); |
||||
|
Assert.AreEqual(factorEvd.Rank, order - 1); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
Assert.AreEqual(1.0, factorEvd.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var resultx = factorEvd.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-10); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixX = factorEvd.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-10); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrixWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new DenseVector(order); |
||||
|
factorEvd.Solve(vectorb, resultx); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-10); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new DenseMatrix(order, order); |
||||
|
factorEvd.Solve(matrixB, matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-10); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,331 @@ |
|||||
|
// <copyright file="UserGramSchmidtTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization |
||||
|
{ |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Double.Factorization; |
||||
|
|
||||
|
public class UserGramSchmidtTests |
||||
|
{ |
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new UserGramSchmidt(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentException] |
||||
|
public void WideMatrixThrowsInvalidMatrixOperationException() |
||||
|
{ |
||||
|
new UserGramSchmidt(new UserDefinedMatrix(3, 4)); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.Q.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.Q.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1,1)] |
||||
|
[Row(2,2)] |
||||
|
[Row(5,5)] |
||||
|
[Row(10,6)] |
||||
|
[Row(50,48)] |
||||
|
[Row(100,98)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
// Make sure the Q has the right dimensions.
|
||||
|
Assert.AreEqual(row, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
// Make sure the R has the right dimensions.
|
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.ColumnCount); |
||||
|
|
||||
|
// Make sure the R factor is upper triangular.
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i > j) |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure the Q*R is the original matrix.
|
||||
|
var matrixQfromR = factorGramSchmidt.Q * factorGramSchmidt.R; |
||||
|
for (var i = 0; i < matrixQfromR.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixQfromR.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixA[i, j], matrixQfromR[i, j], 1.0e-11); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVector(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var resultx = factorGramSchmidt.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixX = factorGramSchmidt.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new UserDefinedVector(order); |
||||
|
factorGramSchmidt.Solve(vectorb,resultx); |
||||
|
|
||||
|
Assert.AreEqual(vectorb.Count, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new UserDefinedMatrix(order, order); |
||||
|
factorGramSchmidt.Solve(matrixB,matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,359 @@ |
|||||
|
// <copyright file="EvdTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization |
||||
|
{ |
||||
|
using System.Numerics; |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using LinearAlgebra.Single; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Single.Factorization; |
||||
|
|
||||
|
public class EvdTests |
||||
|
{ |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new DenseEvd(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(I.RowCount, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorEvd.EValues().Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(Complex.One, factorEvd.EValues()[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A*V = λ*V
|
||||
|
var matrixAv = matrixA * factorEvd.EVectors(); |
||||
|
var matrixLv = factorEvd.EVectors() * factorEvd.D(); |
||||
|
|
||||
|
for (var i = 0; i < matrixAv.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixAv.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixAv[i, j], matrixLv[i, j], 1e-3f); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[Ignore] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.EVectors().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.EVectors().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(order, factorEvd.D().RowCount); |
||||
|
Assert.AreEqual(order, factorEvd.D().ColumnCount); |
||||
|
|
||||
|
// Make sure the A = V*λ*VT
|
||||
|
var matrix = factorEvd.EVectors() * factorEvd.D() * factorEvd.EVectors().Transpose(); |
||||
|
|
||||
|
for (var i = 0; i < matrix.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrix.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrix[i, j], matrixA[i, j], 1e-3f); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankSquare(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Rank, order); |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankOfSquareSingular(int order) |
||||
|
{ |
||||
|
var matrixA = new DenseMatrix(order, order); |
||||
|
matrixA[0, 0] = 1; |
||||
|
matrixA[order - 1, order - 1] = 1; |
||||
|
for (var i = 1; i < order - 1; i++) |
||||
|
{ |
||||
|
matrixA[i, i - 1] = 1; |
||||
|
matrixA[i, i + 1] = 1; |
||||
|
matrixA[i - 1, i] = 1; |
||||
|
matrixA[i + 1, i] = 1; |
||||
|
} |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
Assert.AreEqual(factorEvd.Determinant, 0); |
||||
|
Assert.AreEqual(factorEvd.Rank, order - 1); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorEvd = I.Evd(); |
||||
|
Assert.AreEqual(1.0, factorEvd.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var resultx = factorEvd.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-3f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixX = factorEvd.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1e-2f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorAndSymmetricMatrixWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomDenseVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new DenseVector(order); |
||||
|
factorEvd.Solve(vectorb, resultx); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-3f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorEvd = matrixA.Evd(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new DenseMatrix(order, order); |
||||
|
factorEvd.Solve(matrixB, matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1e-2f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,331 @@ |
|||||
|
// <copyright file="UserGramSchmidtTests.cs" company="Math.NET">
|
||||
|
// 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.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization |
||||
|
{ |
||||
|
using LinearAlgebra.Generic.Factorization; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Single.Factorization; |
||||
|
|
||||
|
public class UserGramSchmidtTests |
||||
|
{ |
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new UserGramSchmidt(null); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentException] |
||||
|
public void WideMatrixThrowsInvalidMatrixOperationException() |
||||
|
{ |
||||
|
new UserGramSchmidt(new UserDefinedMatrix(3, 4)); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
for (var i = 0; i < factorGramSchmidt.Q.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.Q.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.Q[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void IdentityDeterminantIsOne(int order) |
||||
|
{ |
||||
|
var I = UserDefinedMatrix.Identity(order); |
||||
|
var factorGramSchmidt = I.GramSchmidt(); |
||||
|
Assert.AreEqual(1.0, factorGramSchmidt.Determinant); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1,1)] |
||||
|
[Row(2,2)] |
||||
|
[Row(5,5)] |
||||
|
[Row(10,6)] |
||||
|
[Row(50,48)] |
||||
|
[Row(100,98)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
// Make sure the Q has the right dimensions.
|
||||
|
Assert.AreEqual(row, factorGramSchmidt.Q.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.Q.ColumnCount); |
||||
|
|
||||
|
// Make sure the R has the right dimensions.
|
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.RowCount); |
||||
|
Assert.AreEqual(column, factorGramSchmidt.R.ColumnCount); |
||||
|
|
||||
|
// Make sure the R factor is upper triangular.
|
||||
|
for (var i = 0; i < factorGramSchmidt.R.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorGramSchmidt.R.ColumnCount; j++) |
||||
|
{ |
||||
|
if (i > j) |
||||
|
{ |
||||
|
Assert.AreEqual(0.0, factorGramSchmidt.R[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure the Q*R is the original matrix.
|
||||
|
var matrixQfromR = factorGramSchmidt.Q * factorGramSchmidt.R; |
||||
|
for (var i = 0; i < matrixQfromR.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixQfromR.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixA[i, j], matrixQfromR[i, j], 10e-5f); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVector(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var resultx = factorGramSchmidt.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < order; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 10e-5f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrix(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixX = factorGramSchmidt.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 10e-5f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorWhenResultVectorGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new UserDefinedVector(order); |
||||
|
factorGramSchmidt.Solve(vectorb,resultx); |
||||
|
|
||||
|
Assert.AreEqual(vectorb.Count, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 10e-5f); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(4)] |
||||
|
[Row(8)] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorGramSchmidt = matrixA.GramSchmidt(); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new UserDefinedMatrix(order, order); |
||||
|
factorGramSchmidt.Solve(matrixB,matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 10e-5f); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
Loading…
Reference in new issue