mirror of
https://github.com/hanzoai/jube-fork.git
synced 2026-10-30 18:33:33 +00:00
richard-churchman
cc0a2bee21
Included open-source Accord.net library code natively to the solution given the deprecated state of that project. Accord.net will be supported inside Jube for own needs (such as upgrading to .Net 7+). Accord.net libraries modified to compile to .Net 6. Modified Accord.net libraries to better support Newtonsoft Json.net serialisation. Included annotations in the Accord.net libraries to house important dependency statistics. Serialised json can be used independently of Jube to recall a Neural Network. Removal of Binary Serializer totally as a precursor to the projects support of .Net 7+. Extensive bug fixing in crucial functions of Neural Network training and recall. While this part of the system remains in Alpha, it is substantially more fit for purpose.
1099 lines
39 KiB
Plaintext
1099 lines
39 KiB
Plaintext
<#@ template debug="false" hostspecific="true" language="C#" #>
|
|
<#@ output extension=".txt" #>
|
|
<#@ import namespace="Microsoft.VisualStudio.TextTemplating" #>
|
|
<#@ include file="T4Toolbox.tt" #>
|
|
<#
|
|
new DecompositionTemplate("", "Double") .RenderToFile("JaggedEigenvalueDecomposition.cs");
|
|
new DecompositionTemplate("F", "Single") .RenderToFile("JaggedEigenvalueDecompositionF.cs");
|
|
#>
|
|
<#+
|
|
public class DecompositionTemplate : Template
|
|
{
|
|
public string Suffix { get; set; }
|
|
public string T { get; set; }
|
|
|
|
public DecompositionTemplate(string suffix, string dataType) { this.Suffix = suffix; this.T = dataType; }
|
|
|
|
public override string TransformText() {
|
|
#>
|
|
// Accord Math Library
|
|
// The Accord.NET Framework
|
|
// http://accord-framework.net
|
|
//
|
|
// Copyright © César Souza, 2009-2017
|
|
// cesarsouza at gmail.com
|
|
//
|
|
// Original work copyright © Lutz Roeder, 2000
|
|
// Adapted from Mapack for .NET, September 2000
|
|
// Adapted from Mapack for COM and Jama routines
|
|
// http://www.aisto.com/roeder/dotnet
|
|
//
|
|
// This library is free software; you can redistribute it and/or
|
|
// modify it under the terms of the GNU Lesser General Public
|
|
// License as published by the Free Software Foundation; either
|
|
// version 2.1 of the License, or (at your option) any later version.
|
|
//
|
|
// This library is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
// Lesser General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU Lesser General Public
|
|
// License along with this library; if not, write to the Free Software
|
|
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
|
//
|
|
|
|
namespace Accord.Math.Decompositions
|
|
{
|
|
using System;
|
|
using Accord.Math;
|
|
using Accord.Compat;
|
|
|
|
/// <summary>
|
|
/// Determines the eigenvalues and eigenvectors of a real square matrix.
|
|
/// </summary>
|
|
///
|
|
/// <remarks>
|
|
/// <para>
|
|
/// In the mathematical discipline of linear algebra, eigendecomposition
|
|
/// or sometimes spectral decomposition is the factorization of a matrix
|
|
/// into a canonical form, whereby the matrix is represented in terms of
|
|
/// its eigenvalues and eigenvectors.</para>
|
|
/// <para>
|
|
/// If <c>A</c> is symmetric, then <c>A = V * D * V'</c> and <c>A = V * V'</c>
|
|
/// where the eigenvalue matrix <c>D</c> is diagonal and the eigenvector matrix <c>V</c> is orthogonal.
|
|
/// If <c>A</c> is not symmetric, the eigenvalue matrix <c>D</c> is block diagonal
|
|
/// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
|
/// <c>lambda + i*mu</c>, in 2-by-2 blocks, <c>[lambda, mu; -mu, lambda]</c>.
|
|
/// The columns of <c>V</c> represent the eigenvectors in the sense that <c>A * V = V * D</c>.
|
|
/// The matrix V may be badly conditioned, or even singular, so the validity of the equation
|
|
/// <c>A = V * D * inverse(V)</c> depends upon the condition of <c>V</c>.
|
|
/// </para>
|
|
/// </remarks>
|
|
///
|
|
public sealed class JaggedEigenvalueDecomposition<#=Suffix#> : ICloneable
|
|
{
|
|
private int n; // matrix dimension
|
|
private <#=T#>[] d, e; // storage of eigenvalues.
|
|
private <#=T#>[][] V; // storage of eigenvectors.
|
|
private <#=T#>[][] H; // storage of nonsymmetric Hessenberg form.
|
|
private <#=T#>[] ort; // storage for nonsymmetric algorithm.
|
|
private bool symmetric;
|
|
|
|
private int? rank;
|
|
private <#=T#>[][] diagonalMatrix;
|
|
|
|
private const <#=T#> eps = 2 * Constants.<#=T#>Epsilon;
|
|
|
|
|
|
/// <summary>
|
|
/// Returns the effective numerical matrix rank.
|
|
/// </summary>
|
|
///
|
|
/// <value>Number of non-negligible eigen values.</value>
|
|
///
|
|
public int Rank
|
|
{
|
|
get
|
|
{
|
|
if (this.rank.HasValue)
|
|
return this.rank.Value;
|
|
|
|
<#=T#> tol = n * d[0] * eps;
|
|
|
|
int r = 0;
|
|
for (int i = 0; i < d.Length; i++)
|
|
if (d[i] > tol) r++;
|
|
|
|
return (int)(this.rank = r);
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Construct an eigenvalue decomposition.</summary>
|
|
/// <param name="value">
|
|
/// The matrix to be decomposed.</param>
|
|
/// <param name="inPlace">
|
|
/// Pass <see langword="true"/> to perform the decomposition in place. The matrix
|
|
/// <paramref name="value"/> will be destroyed in the process, resulting in less
|
|
/// memory comsumption.</param>
|
|
/// <param name="sort">
|
|
/// Pass <see langword="true"/> to sort the eigenvalues and eigenvectors at the end
|
|
/// of the decomposition.</param>
|
|
///
|
|
public JaggedEigenvalueDecomposition<#=Suffix#>(<#=T#>[][] value, bool inPlace = false, bool sort = false)
|
|
: this(value, assumeSymmetric: value.IsSymmetric(), inPlace: inPlace, sort: sort)
|
|
{
|
|
}
|
|
|
|
/// <summary>
|
|
/// Construct an eigenvalue decomposition.</summary>
|
|
/// <param name="value">
|
|
/// The matrix to be decomposed.</param>
|
|
/// <param name="assumeSymmetric">
|
|
/// Defines if the matrix should be assumed as being symmetric
|
|
/// regardless if it is or not. Default is <see langword="false"/>.</param>
|
|
/// <param name="inPlace">
|
|
/// Pass <see langword="true"/> to perform the decomposition in place. The matrix
|
|
/// <paramref name="value"/> will be destroyed in the process, resulting in less
|
|
/// memory comsumption.</param>
|
|
/// <param name="sort">
|
|
/// Pass <see langword="true"/> to sort the eigenvalues and eigenvectors at the end
|
|
/// of the decomposition.</param>
|
|
///
|
|
public JaggedEigenvalueDecomposition<#=Suffix#>(<#=T#>[][] value, bool assumeSymmetric,
|
|
bool inPlace = false, bool sort = false)
|
|
{
|
|
if (value == null)
|
|
throw new ArgumentNullException("value", "Matrix cannot be null.");
|
|
|
|
if (value.Length != value[0].Length)
|
|
throw new ArgumentException("Matrix is not a square matrix.", "value");
|
|
|
|
n = value.Length;
|
|
V = new <#=T#>[n][];
|
|
for (int i = 0; i < n; i++)
|
|
V[i] = new <#=T#>[n];
|
|
d = new <#=T#>[n];
|
|
e = new <#=T#>[n];
|
|
|
|
|
|
this.symmetric = assumeSymmetric;
|
|
|
|
if (this.symmetric)
|
|
{
|
|
V = inPlace ? value : (<#=T#>[][])value.MemberwiseClone();
|
|
|
|
// Tridiagonalize.
|
|
this.tred2();
|
|
|
|
// Diagonalize.
|
|
this.tql2();
|
|
}
|
|
else
|
|
{
|
|
H = inPlace ? value : (<#=T#>[][])value.MemberwiseClone();
|
|
|
|
ort = new <#=T#>[n];
|
|
|
|
// Reduce to Hessenberg form.
|
|
this.orthes();
|
|
|
|
// Reduce Hessenberg to real Schur form.
|
|
this.hqr2();
|
|
}
|
|
|
|
if (sort)
|
|
{
|
|
// Sort eigenvalues and vectors in descending order
|
|
var idx = Vector.Range(n);
|
|
Array.Sort(idx, (i, j) =>
|
|
{
|
|
if (Math.Abs(d[i]) == Math.Abs(d[j]))
|
|
return -Math.Abs(e[i]).CompareTo(Math.Abs(e[j]));
|
|
return -Math.Abs(d[i]).CompareTo(Math.Abs(d[j]));
|
|
});
|
|
|
|
this.d = this.d.Get(idx);
|
|
this.e = this.e.Get(idx);
|
|
this.V = this.V.Get(null, idx);
|
|
}
|
|
}
|
|
|
|
|
|
/// <summary>Returns the real parts of the eigenvalues.</summary>
|
|
public <#=T#>[] RealEigenvalues
|
|
{
|
|
get { return this.d; }
|
|
}
|
|
|
|
/// <summary>Returns the imaginary parts of the eigenvalues.</summary>
|
|
public <#=T#>[] ImaginaryEigenvalues
|
|
{
|
|
get { return this.e; }
|
|
}
|
|
|
|
/// <summary>Returns the eigenvector matrix.</summary>
|
|
public <#=T#>[][] Eigenvectors
|
|
{
|
|
get { return this.V; }
|
|
}
|
|
|
|
/// <summary>Returns the block diagonal eigenvalue matrix.</summary>
|
|
public <#=T#>[][] DiagonalMatrix
|
|
{
|
|
get
|
|
{
|
|
if (this.diagonalMatrix != null)
|
|
return this.diagonalMatrix;
|
|
|
|
var x = new <#=T#>[n][];
|
|
for (int i = 0; i < n; i++)
|
|
x[i] = new <#=T#>[n];
|
|
|
|
for (int i = 0; i < n; i++)
|
|
{
|
|
for (int j = 0; j < n; j++)
|
|
x[i][j] = 0;
|
|
|
|
x[i][i] = d[i];
|
|
if (e[i] > 0)
|
|
{
|
|
x[i][i + 1] = e[i];
|
|
}
|
|
else if (e[i] < 0)
|
|
{
|
|
x[i][i - 1] = e[i];
|
|
}
|
|
}
|
|
|
|
return this.diagonalMatrix = x;
|
|
}
|
|
}
|
|
|
|
|
|
#region Private methods
|
|
private void tred2()
|
|
{
|
|
// Symmetric Householder reduction to tridiagonal form.
|
|
// This is derived from the Algol procedures tred2 by Bowdler, Martin, Reinsch, and Wilkinson,
|
|
// Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding Fortran subroutine in EISPACK.
|
|
for (int j = 0; j < n; j++)
|
|
d[j] = V[n - 1][j];
|
|
|
|
// Householder reduction to tridiagonal form.
|
|
for (int i = n - 1; i > 0; i--)
|
|
{
|
|
// Scale to avoid under/overflow.
|
|
<#=T#> scale = 0;
|
|
<#=T#> h = 0;
|
|
for (int k = 0; k < i; k++)
|
|
scale = scale + System.Math.Abs(d[k]);
|
|
|
|
if (scale == 0)
|
|
{
|
|
e[i] = d[i - 1];
|
|
for (int j = 0; j < i; j++)
|
|
{
|
|
d[j] = V[i - 1][j];
|
|
V[i][j] = 0;
|
|
V[j][i] = 0;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Generate Householder vector.
|
|
for (int k = 0; k < i; k++)
|
|
{
|
|
d[k] /= scale;
|
|
h += d[k] * d[k];
|
|
}
|
|
|
|
<#=T#> f = d[i - 1];
|
|
<#=T#> g = (<#=T#>)System.Math.Sqrt(h);
|
|
if (f > 0) g = -g;
|
|
|
|
e[i] = scale * g;
|
|
h = h - f * g;
|
|
d[i - 1] = f - g;
|
|
for (int j = 0; j < i; j++)
|
|
e[j] = 0;
|
|
|
|
// Apply similarity transformation to remaining columns.
|
|
for (int j = 0; j < i; j++)
|
|
{
|
|
f = d[j];
|
|
V[j][i] = f;
|
|
g = e[j] + V[j][j] * f;
|
|
for (int k = j + 1; k <= i - 1; k++)
|
|
{
|
|
g += V[k][j] * d[k];
|
|
e[k] += V[k][j] * f;
|
|
}
|
|
e[j] = g;
|
|
}
|
|
|
|
f = 0;
|
|
for (int j = 0; j < i; j++)
|
|
{
|
|
e[j] /= h;
|
|
f += e[j] * d[j];
|
|
}
|
|
|
|
<#=T#> hh = f / (h + h);
|
|
for (int j = 0; j < i; j++)
|
|
e[j] -= hh * d[j];
|
|
|
|
for (int j = 0; j < i; j++)
|
|
{
|
|
f = d[j];
|
|
g = e[j];
|
|
for (int k = j; k <= i - 1; k++)
|
|
V[k][j] -= (f * e[k] + g * d[k]);
|
|
|
|
d[j] = V[i - 1][j];
|
|
V[i][j] = 0;
|
|
}
|
|
}
|
|
d[i] = h;
|
|
}
|
|
|
|
// Accumulate transformations.
|
|
for (int i = 0; i < n - 1; i++)
|
|
{
|
|
V[n - 1][i] = V[i][i];
|
|
V[i][i] = 1;
|
|
<#=T#> h = d[i + 1];
|
|
if (h != 0)
|
|
{
|
|
for (int k = 0; k <= i; k++)
|
|
d[k] = V[k][i + 1] / h;
|
|
|
|
for (int j = 0; j <= i; j++)
|
|
{
|
|
<#=T#> g = 0;
|
|
for (int k = 0; k <= i; k++)
|
|
g += V[k][i + 1] * V[k][j];
|
|
for (int k = 0; k <= i; k++)
|
|
V[k][j] -= g * d[k];
|
|
}
|
|
}
|
|
|
|
for (int k = 0; k <= i; k++)
|
|
V[k][i + 1] = 0;
|
|
}
|
|
|
|
for (int j = 0; j < n; j++)
|
|
{
|
|
d[j] = V[n - 1][j];
|
|
V[n - 1][j] = 0;
|
|
}
|
|
|
|
V[n - 1][n - 1] = 1;
|
|
e[0] = 0;
|
|
}
|
|
|
|
private void tql2()
|
|
{
|
|
// Symmetric tridiagonal QL algorithm.
|
|
// 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.
|
|
for (int i = 1; i < n; i++)
|
|
e[i - 1] = e[i];
|
|
|
|
e[n - 1] = 0;
|
|
|
|
<#=T#> f = 0;
|
|
<#=T#> tst1 = 0;
|
|
<#=T#> eps = 2 * Constants.<#=T#>Epsilon;
|
|
|
|
for (int l = 0; l < n; l++)
|
|
{
|
|
// Find small subdiagonal element.
|
|
tst1 = System.Math.Max(tst1, System.Math.Abs(d[l]) + System.Math.Abs(e[l]));
|
|
int m = l;
|
|
while (m < n)
|
|
{
|
|
if (System.Math.Abs(e[m]) <= eps * tst1)
|
|
break;
|
|
m++;
|
|
}
|
|
|
|
// If m == l, d[l] is an eigenvalue, otherwise, iterate.
|
|
if (m > l)
|
|
{
|
|
int iter = 0;
|
|
do
|
|
{
|
|
iter = iter + 1; // (Could check iteration count here.)
|
|
|
|
// Compute implicit shift
|
|
<#=T#> g = d[l];
|
|
<#=T#> p = (d[l + 1] - g) / (2 * e[l]);
|
|
<#=T#> r = Accord.Math.Tools.Hypotenuse(p, 1);
|
|
if (p < 0)
|
|
{
|
|
r = -r;
|
|
}
|
|
|
|
d[l] = e[l] / (p + r);
|
|
d[l + 1] = e[l] * (p + r);
|
|
<#=T#> dl1 = d[l + 1];
|
|
<#=T#> h = g - d[l];
|
|
for (int i = l + 2; i < n; i++)
|
|
{
|
|
d[i] -= h;
|
|
}
|
|
|
|
f = f + h;
|
|
|
|
// Implicit QL transformation.
|
|
p = d[m];
|
|
<#=T#> c = 1;
|
|
<#=T#> c2 = c;
|
|
<#=T#> c3 = c;
|
|
<#=T#> el1 = e[l + 1];
|
|
<#=T#> s = 0;
|
|
<#=T#> s2 = 0;
|
|
for (int i = m - 1; i >= l; i--)
|
|
{
|
|
c3 = c2;
|
|
c2 = c;
|
|
s2 = s;
|
|
g = c * e[i];
|
|
h = c * p;
|
|
r = Accord.Math.Tools.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 (int k = 0; k < n; k++)
|
|
{
|
|
h = V[k][i + 1];
|
|
V[k][i + 1] = s * V[k][i] + c * h;
|
|
V[k][i] = c * V[k][i] - s * h;
|
|
}
|
|
}
|
|
|
|
p = -s * s2 * c3 * el1 * e[l] / dl1;
|
|
e[l] = s * p;
|
|
d[l] = c * p;
|
|
|
|
// Check for convergence.
|
|
}
|
|
while (System.Math.Abs(e[l]) > eps * tst1);
|
|
}
|
|
d[l] = d[l] + f;
|
|
e[l] = 0;
|
|
}
|
|
|
|
// Sort eigenvalues and corresponding vectors.
|
|
for (int i = 0; i < n - 1; i++)
|
|
{
|
|
int k = i;
|
|
<#=T#> p = d[i];
|
|
for (int j = i + 1; j < n; j++)
|
|
{
|
|
if (d[j] < p)
|
|
{
|
|
k = j;
|
|
p = d[j];
|
|
}
|
|
}
|
|
|
|
if (k != i)
|
|
{
|
|
d[k] = d[i];
|
|
d[i] = p;
|
|
for (int j = 0; j < n; j++)
|
|
{
|
|
p = V[j][i];
|
|
V[j][i] = V[j][k];
|
|
V[j][k] = p;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
private void orthes()
|
|
{
|
|
// Nonsymmetric reduction to Hessenberg form.
|
|
// 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.
|
|
int low = 0;
|
|
int high = n - 1;
|
|
|
|
for (int m = low + 1; m <= high - 1; m++)
|
|
{
|
|
// Scale column.
|
|
|
|
<#=T#> scale = 0;
|
|
for (int i = m; i <= high; i++)
|
|
scale = scale + System.Math.Abs(H[i][m - 1]);
|
|
|
|
if (scale != 0)
|
|
{
|
|
// Compute Householder transformation.
|
|
<#=T#> h = 0;
|
|
for (int i = high; i >= m; i--)
|
|
{
|
|
ort[i] = H[i][m - 1] / scale;
|
|
h += ort[i] * ort[i];
|
|
}
|
|
|
|
<#=T#> g = (<#=T#>)System.Math.Sqrt(h);
|
|
if (ort[m] > 0) g = -g;
|
|
|
|
h = h - ort[m] * g;
|
|
ort[m] = ort[m] - g;
|
|
|
|
// Apply Householder similarity transformation
|
|
// H = (I - u * u' / h) * H * (I - u * u') / h)
|
|
for (int j = m; j < n; j++)
|
|
{
|
|
<#=T#> f = 0;
|
|
for (int i = high; i >= m; i--)
|
|
f += ort[i] * H[i][j];
|
|
|
|
f = f / h;
|
|
for (int i = m; i <= high; i++)
|
|
H[i][j] -= f * ort[i];
|
|
}
|
|
|
|
for (int i = 0; i <= high; i++)
|
|
{
|
|
<#=T#> f = 0;
|
|
for (int j = high; j >= m; j--)
|
|
f += ort[j] * H[i][j];
|
|
|
|
f = f / h;
|
|
for (int j = m; j <= high; j++)
|
|
H[i][j] -= f * ort[j];
|
|
}
|
|
|
|
ort[m] = scale * ort[m];
|
|
H[m][m - 1] = scale * g;
|
|
}
|
|
}
|
|
|
|
// Accumulate transformations (Algol's ortran).
|
|
for (int i = 0; i < n; i++)
|
|
for (int j = 0; j < n; j++)
|
|
V[i][j] = (i == j ? 1 : 0);
|
|
|
|
for (int m = high - 1; m >= low + 1; m--)
|
|
{
|
|
if (H[m][m - 1] != 0)
|
|
{
|
|
for (int i = m + 1; i <= high; i++)
|
|
ort[i] = H[i][m - 1];
|
|
|
|
for (int j = m; j <= high; j++)
|
|
{
|
|
<#=T#> g = 0;
|
|
for (int i = m; i <= high; i++)
|
|
g += ort[i] * V[i][j];
|
|
|
|
// Double division avoids possible underflow.
|
|
g = (g / ort[m]) / H[m][m - 1];
|
|
for (int i = m; i <= high; i++)
|
|
V[i][j] += g * ort[i];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
private static void cdiv(<#=T#> xr, <#=T#> xi, <#=T#> yr, <#=T#> yi,
|
|
out <#=T#> cdivr, out <#=T#> cdivi)
|
|
{
|
|
// Complex scalar division.
|
|
<#=T#> r;
|
|
<#=T#> d;
|
|
if (System.Math.Abs(yr) > System.Math.Abs(yi))
|
|
{
|
|
r = yi / yr;
|
|
d = yr + r * yi;
|
|
cdivr = (xr + r * xi) / d;
|
|
cdivi = (xi - r * xr) / d;
|
|
}
|
|
else
|
|
{
|
|
r = yr / yi;
|
|
d = yi + r * yr;
|
|
cdivr = (r * xr + xi) / d;
|
|
cdivi = (r * xi - xr) / d;
|
|
}
|
|
}
|
|
|
|
private void hqr2()
|
|
{
|
|
// Nonsymmetric reduction from Hessenberg to real Schur form.
|
|
// 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.
|
|
int nn = this.n;
|
|
int n = nn - 1;
|
|
int low = 0;
|
|
int high = nn - 1;
|
|
<#=T#> eps = 2 * Constants.<#=T#>Epsilon;
|
|
<#=T#> exshift = 0;
|
|
<#=T#> p = 0;
|
|
<#=T#> q = 0;
|
|
<#=T#> r = 0;
|
|
<#=T#> s = 0;
|
|
<#=T#> z = 0;
|
|
<#=T#> t;
|
|
<#=T#> w;
|
|
<#=T#> x;
|
|
<#=T#> y;
|
|
|
|
// Store roots isolated by balanc and compute matrix norm
|
|
<#=T#> norm = 0;
|
|
for (int i = 0; i < nn; i++)
|
|
{
|
|
if (i < low | i > high)
|
|
{
|
|
d[i] = H[i][i];
|
|
e[i] = 0;
|
|
}
|
|
|
|
for (int j = System.Math.Max(i - 1, 0); j < nn; j++)
|
|
norm = norm + System.Math.Abs(H[i][j]);
|
|
}
|
|
|
|
// Outer loop over eigenvalue index
|
|
int iter = 0;
|
|
while (n >= low)
|
|
{
|
|
// Look for single small sub-diagonal element
|
|
int l = n;
|
|
while (l > low)
|
|
{
|
|
s = System.Math.Abs(H[l - 1][l - 1]) + System.Math.Abs(H[l][l]);
|
|
|
|
if (s == 0)
|
|
s = norm;
|
|
|
|
if (<#=T#>.IsNaN(s))
|
|
break;
|
|
|
|
if (System.Math.Abs(H[l][l - 1]) < eps * s)
|
|
break;
|
|
|
|
l--;
|
|
}
|
|
|
|
// Check for convergence
|
|
if (l == n)
|
|
{
|
|
// One root found
|
|
H[n][n] = H[n][n] + exshift;
|
|
d[n] = H[n][n];
|
|
e[n] = 0;
|
|
n--;
|
|
iter = 0;
|
|
}
|
|
else if (l == n - 1)
|
|
{
|
|
// Two roots found
|
|
w = H[n][n - 1] * H[n - 1][n];
|
|
p = (H[n - 1][n - 1] - H[n][n]) / 2;
|
|
q = p * p + w;
|
|
z = (<#=T#>)System.Math.Sqrt(System.Math.Abs(q));
|
|
H[n][n] = H[n][n] + exshift;
|
|
H[n - 1][n - 1] = H[n - 1][n - 1] + exshift;
|
|
x = H[n][n];
|
|
|
|
if (q >= 0)
|
|
{
|
|
// Real pair
|
|
z = (p >= 0) ? (p + z) : (p - z);
|
|
d[n - 1] = x + z;
|
|
d[n] = d[n - 1];
|
|
if (z != 0)
|
|
d[n] = x - w / z;
|
|
e[n - 1] = 0;
|
|
e[n] = 0;
|
|
x = H[n][n - 1];
|
|
s = System.Math.Abs(x) + System.Math.Abs(z);
|
|
p = x / s;
|
|
q = z / s;
|
|
r = (<#=T#>)System.Math.Sqrt(p * p + q * q);
|
|
p = p / r;
|
|
q = q / r;
|
|
|
|
// Row modification
|
|
for (int j = n - 1; j < nn; j++)
|
|
{
|
|
z = H[n - 1][j];
|
|
H[n - 1][j] = q * z + p * H[n][j];
|
|
H[n][j] = q * H[n][j] - p * z;
|
|
}
|
|
|
|
// Column modification
|
|
for (int i = 0; i <= n; i++)
|
|
{
|
|
z = H[i][n - 1];
|
|
H[i][n - 1] = q * z + p * H[i][n];
|
|
H[i][n] = q * H[i][n] - p * z;
|
|
}
|
|
|
|
// Accumulate transformations
|
|
for (int i = low; i <= high; i++)
|
|
{
|
|
z = V[i][n - 1];
|
|
V[i][n - 1] = q * z + p * V[i][n];
|
|
V[i][n] = q * V[i][n] - p * z;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Complex pair
|
|
d[n - 1] = x + p;
|
|
d[n] = x + p;
|
|
e[n - 1] = z;
|
|
e[n] = -z;
|
|
}
|
|
|
|
n = n - 2;
|
|
iter = 0;
|
|
}
|
|
else
|
|
{
|
|
// No convergence yet
|
|
|
|
// Form shift
|
|
x = H[n][n];
|
|
y = 0;
|
|
w = 0;
|
|
if (l < n)
|
|
{
|
|
y = H[n - 1][n - 1];
|
|
w = H[n][n - 1] * H[n - 1][n];
|
|
}
|
|
|
|
// Wilkinson's original ad hoc shift
|
|
if (iter == 10)
|
|
{
|
|
exshift += x;
|
|
for (int i = low; i <= n; i++)
|
|
H[i][i] -= x;
|
|
|
|
s = System.Math.Abs(H[n][n - 1]) + System.Math.Abs(H[n - 1][n - 2]);
|
|
x = y = (<#=T#>)0.75 * s;
|
|
w = (<#=T#>)(-0.4375) * s * s;
|
|
}
|
|
|
|
// MATLAB's new ad hoc shift
|
|
if (iter == 30)
|
|
{
|
|
s = (y - x) / 2;
|
|
s = s * s + w;
|
|
if (s > 0)
|
|
{
|
|
s = (<#=T#>)System.Math.Sqrt(s);
|
|
if (y < x) s = -s;
|
|
s = x - w / ((y - x) / 2 + s);
|
|
for (int i = low; i <= n; i++)
|
|
H[i][i] -= s;
|
|
exshift += s;
|
|
x = y = w = (<#=T#>)0.964;
|
|
}
|
|
}
|
|
|
|
iter = iter + 1;
|
|
|
|
// Look for two consecutive small sub-diagonal elements
|
|
int m = n - 2;
|
|
while (m >= l)
|
|
{
|
|
z = H[m][m];
|
|
r = x - z;
|
|
s = y - z;
|
|
p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
|
|
q = H[m + 1][m + 1] - z - r - s;
|
|
r = H[m + 2][m + 1];
|
|
s = System.Math.Abs(p) + System.Math.Abs(q) + System.Math.Abs(r);
|
|
p = p / s;
|
|
q = q / s;
|
|
r = r / s;
|
|
if (m == l)
|
|
break;
|
|
if (System.Math.Abs(H[m][m - 1]) * (System.Math.Abs(q) + System.Math.Abs(r)) < eps * (System.Math.Abs(p) * (System.Math.Abs(H[m - 1][m - 1]) + System.Math.Abs(z) + System.Math.Abs(H[m + 1][m + 1]))))
|
|
break;
|
|
m--;
|
|
}
|
|
|
|
for (int i = m + 2; i <= n; i++)
|
|
{
|
|
H[i][i - 2] = 0;
|
|
if (i > m + 2)
|
|
H[i][i - 3] = 0;
|
|
}
|
|
|
|
// Double QR step involving rows l:n and columns m:n
|
|
for (int k = m; k <= n - 1; k++)
|
|
{
|
|
bool notlast = (k != n - 1);
|
|
if (k != m)
|
|
{
|
|
p = H[k][k - 1];
|
|
q = H[k + 1][k - 1];
|
|
r = (notlast ? H[k + 2][k - 1] : 0);
|
|
x = System.Math.Abs(p) + System.Math.Abs(q) + System.Math.Abs(r);
|
|
if (x != 0)
|
|
{
|
|
p = p / x;
|
|
q = q / x;
|
|
r = r / x;
|
|
}
|
|
}
|
|
|
|
if (x == 0)
|
|
break;
|
|
|
|
s = (<#=T#>)System.Math.Sqrt(p * p + q * q + r * r);
|
|
if (p < 0) s = -s;
|
|
|
|
if (s != 0)
|
|
{
|
|
if (k != m)
|
|
H[k][k - 1] = -s * x;
|
|
else
|
|
if (l != m)
|
|
H[k][k - 1] = -H[k][k - 1];
|
|
|
|
p = p + s;
|
|
x = p / s;
|
|
y = q / s;
|
|
z = r / s;
|
|
q = q / p;
|
|
r = r / p;
|
|
|
|
// Row modification
|
|
for (int j = k; j < nn; j++)
|
|
{
|
|
p = H[k][j] + q * H[k + 1][j];
|
|
if (notlast)
|
|
{
|
|
p = p + r * H[k + 2][j];
|
|
H[k + 2][j] = H[k + 2][j] - p * z;
|
|
}
|
|
|
|
H[k][j] = H[k][j] - p * x;
|
|
H[k + 1][j] = H[k + 1][j] - p * y;
|
|
}
|
|
|
|
// Column modification
|
|
for (int i = 0; i <= System.Math.Min(n, k + 3); i++)
|
|
{
|
|
p = x * H[i][k] + y * H[i][k + 1];
|
|
if (notlast)
|
|
{
|
|
p = p + z * H[i][k + 2];
|
|
H[i][k + 2] = H[i][k + 2] - p * r;
|
|
}
|
|
|
|
H[i][k] = H[i][k] - p;
|
|
H[i][k + 1] = H[i][k + 1] - p * q;
|
|
}
|
|
|
|
// Accumulate transformations
|
|
for (int i = low; i <= high; i++)
|
|
{
|
|
p = x * V[i][k] + y * V[i][k + 1];
|
|
if (notlast)
|
|
{
|
|
p = p + z * V[i][k + 2];
|
|
V[i][k + 2] = V[i][k + 2] - p * r;
|
|
}
|
|
|
|
V[i][k] = V[i][k] - p;
|
|
V[i][k + 1] = V[i][k + 1] - p * q;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Backsubstitute to find vectors of upper triangular form
|
|
if (norm == 0)
|
|
{
|
|
return;
|
|
}
|
|
|
|
for (n = nn - 1; n >= 0; n--)
|
|
{
|
|
p = d[n];
|
|
q = e[n];
|
|
|
|
// Real vector
|
|
if (q == 0)
|
|
{
|
|
int l = n;
|
|
H[n][n] = 1;
|
|
for (int i = n - 1; i >= 0; i--)
|
|
{
|
|
w = H[i][i] - p;
|
|
r = 0;
|
|
for (int j = l; j <= n; j++)
|
|
r = r + H[i][j] * H[j][n];
|
|
|
|
if (e[i] < 0)
|
|
{
|
|
z = w;
|
|
s = r;
|
|
}
|
|
else
|
|
{
|
|
l = i;
|
|
if (e[i] == 0)
|
|
{
|
|
H[i][n] = (w != 0) ? (-r / w) : (-r / (eps * norm));
|
|
}
|
|
else
|
|
{
|
|
// Solve real equations
|
|
x = H[i][i + 1];
|
|
y = H[i + 1][i];
|
|
q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
|
|
t = (x * s - z * r) / q;
|
|
H[i][n] = t;
|
|
H[i + 1][n] = (System.Math.Abs(x) > System.Math.Abs(z)) ? ((-r - w * t) / x) : ((-s - y * t) / z);
|
|
}
|
|
|
|
// Overflow control
|
|
t = System.Math.Abs(H[i][n]);
|
|
if ((eps * t) * t > 1)
|
|
for (int j = i; j <= n; j++)
|
|
H[j][n] = H[j][n] / t;
|
|
}
|
|
}
|
|
}
|
|
else if (q < 0)
|
|
{
|
|
// Complex vector
|
|
int l = n - 1;
|
|
|
|
// Last vector component imaginary so matrix is triangular
|
|
if (System.Math.Abs(H[n][n - 1]) > System.Math.Abs(H[n - 1][n]))
|
|
{
|
|
H[n - 1][n - 1] = q / H[n][n - 1];
|
|
H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1];
|
|
}
|
|
else
|
|
{
|
|
cdiv(0, -H[n - 1][n], H[n - 1][n - 1] - p, q, out H[n - 1][n - 1], out H[n - 1][n]);
|
|
}
|
|
|
|
H[n][n - 1] = 0;
|
|
H[n][n] = 1;
|
|
for (int i = n - 2; i >= 0; i--)
|
|
{
|
|
<#=T#> ra, sa, vr, vi;
|
|
ra = 0;
|
|
sa = 0;
|
|
for (int j = l; j <= n; j++)
|
|
{
|
|
ra = ra + H[i][j] * H[j][n - 1];
|
|
sa = sa + H[i][j] * H[j][n];
|
|
}
|
|
|
|
w = H[i][i] - p;
|
|
|
|
if (e[i] < 0)
|
|
{
|
|
z = w;
|
|
r = ra;
|
|
s = sa;
|
|
}
|
|
else
|
|
{
|
|
l = i;
|
|
if (e[i] == 0)
|
|
{
|
|
cdiv(-ra, -sa, w, q, out H[i][n - 1], out H[i][n]);
|
|
}
|
|
else
|
|
{
|
|
// Solve complex equations
|
|
x = H[i][i + 1];
|
|
y = H[i + 1][i];
|
|
vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
|
|
vi = (d[i] - p) * 2 * q;
|
|
if (vr == 0 & vi == 0)
|
|
vr = eps * norm * (System.Math.Abs(w) + System.Math.Abs(q) + System.Math.Abs(x) + System.Math.Abs(y) + System.Math.Abs(z));
|
|
cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi, out H[i][n - 1], out H[i][n]);
|
|
if (System.Math.Abs(x) > (System.Math.Abs(z) + System.Math.Abs(q)))
|
|
{
|
|
H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x;
|
|
H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x;
|
|
}
|
|
else
|
|
{
|
|
cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q, out H[i + 1][n - 1], out H[i + 1][n]);
|
|
}
|
|
}
|
|
|
|
// Overflow control
|
|
t = System.Math.Max(System.Math.Abs(H[i][n - 1]), System.Math.Abs(H[i][n]));
|
|
if ((eps * t) * t > 1)
|
|
{
|
|
for (int j = i; j <= n; j++)
|
|
{
|
|
H[j][n - 1] = H[j][n - 1] / t;
|
|
H[j][n] = H[j][n] / t;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Vectors of isolated roots
|
|
for (int i = 0; i < nn; i++)
|
|
if (i < low | i > high)
|
|
for (int j = i; j < nn; j++)
|
|
V[i][j] = H[i][j];
|
|
|
|
// Back transformation to get eigenvectors of original matrix
|
|
for (int j = nn - 1; j >= low; j--)
|
|
{
|
|
for (int i = low; i <= high; i++)
|
|
{
|
|
z = 0;
|
|
for (int k = low; k <= System.Math.Min(j, high); k++)
|
|
z = z + V[i][k] * H[k][j];
|
|
V[i][j] = z;
|
|
}
|
|
}
|
|
}
|
|
#endregion
|
|
|
|
/// <summary>
|
|
/// Reverses the decomposition, reconstructing the original matrix <c>X</c>.
|
|
/// </summary>
|
|
///
|
|
public <#=T#>[][] Reverse()
|
|
{
|
|
return V.DotWithDiagonal(d).Divide(V);
|
|
}
|
|
|
|
|
|
#region ICloneable Members
|
|
|
|
private JaggedEigenvalueDecomposition<#=Suffix#>()
|
|
{
|
|
}
|
|
|
|
/// <summary>
|
|
/// Creates a new object that is a copy of the current instance.
|
|
/// </summary>
|
|
/// <returns>
|
|
/// A new object that is a copy of this instance.
|
|
/// </returns>
|
|
public object Clone()
|
|
{
|
|
var clone = new JaggedEigenvalueDecomposition<#=Suffix#>();
|
|
clone.d = (<#=T#>[])d.Clone();
|
|
clone.e = (<#=T#>[])e.Clone();
|
|
clone.H = H.MemberwiseClone();
|
|
clone.n = n;
|
|
clone.ort = (<#=T#>[])ort;
|
|
clone.symmetric = symmetric;
|
|
clone.V = V.MemberwiseClone();
|
|
return clone;
|
|
}
|
|
|
|
#endregion
|
|
|
|
}
|
|
}
|
|
|
|
<#+
|
|
return GenerationEnvironment.ToString();
|
|
}
|
|
}
|
|
#> |