diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs index b873c9fc..bb65e2bb 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -98,724 +98,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { } - /// - /// Reduces a complex Hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations. - /// - /// Source matrix to reduce - /// Output: Arrays for internal storage of real parts of eigenvalues - /// Output: Arrays for internal storage of imaginary parts of eigenvalues - /// Output: Arrays that contains further information about the transformations. - /// Order of initial matrix - /// 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. - internal 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*order + 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[k*order + i].Real) + Math.Abs(matrixA[k*order + i].Imaginary); - } - - if (scale == 0.0) - { - tau[i - 1] = Complex.One; - e[i] = 0.0; - } - else - { - for (var k = 0; k < i; k++) - { - matrixA[k*order + i] /= scale; - h += matrixA[k*order + i].MagnitudeSquared(); - } - - Complex g = Math.Sqrt(h); - e[i] = scale*g.Real; - - Complex temp; - var im1Oi = (i - 1)*order + i; - var f = matrixA[im1Oi]; - if (f.Magnitude != 0) - { - temp = -(matrixA[im1Oi].Conjugate()*tau[i].Conjugate())/f.Magnitude; - h += f.Magnitude*g.Real; - g = 1.0 + (g/f.Magnitude); - matrixA[im1Oi] *= g; - } - else - { - temp = -tau[i].Conjugate(); - matrixA[im1Oi] = g; - } - - if ((f.Magnitude == 0) || (i != 1)) - { - f = Complex.Zero; - for (var j = 0; j < i; j++) - { - var tmp = Complex.Zero; - var jO = j*order; - // Form element of A*U. - for (var k = 0; k <= j; k++) - { - tmp += matrixA[k*order + j]*matrixA[k*order + i].Conjugate(); - } - - for (var k = j + 1; k <= i - 1; k++) - { - tmp += matrixA[jO + k].Conjugate()*matrixA[k*order + i].Conjugate(); - } - - // Form element of P - tau[j] = tmp/h; - f += (tmp/h)*matrixA[jO + i]; - } - - hh = f.Real/(h + h); - - // Form the reduced A. - for (var j = 0; j < i; j++) - { - f = matrixA[j*order + i].Conjugate(); - g = tau[j] - (hh*f); - tau[j] = g.Conjugate(); - - for (var k = 0; k <= j; k++) - { - matrixA[k*order + j] -= (f*tau[k]) + (g*matrixA[k*order + i]); - } - } - } - - for (var k = 0; k < i; k++) - { - matrixA[k*order + i] *= scale; - } - - tau[i - 1] = temp.Conjugate(); - } - - hh = d[i]; - d[i] = matrixA[i*order + i].Real; - matrixA[i*order + i] = new Complex(hh, scale*Math.Sqrt(h)); - } - - hh = d[0]; - d[0] = matrixA[0].Real; - matrixA[0] = hh; - e[0] = 0.0; - } - - /// - /// Symmetric tridiagonal QL algorithm. - /// - /// Data array of matrix V (eigenvectors) - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - /// - internal 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.DoublePrecision; - 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 NonConvergenceException(); - } - } 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; - } - } - } - } - - /// - /// Determines eigenvectors by undoing the symmetric tridiagonalize transformation - /// - /// Data array of matrix V (eigenvectors) - /// Previously tridiagonalized matrix by . - /// Contains further information about the transformations - /// Input matrix order - /// 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. - internal 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*order + 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[k*order + i]; - } - - s = (s/h)/h; - - for (var k = 0; k < i; k++) - { - dataEv[(j*order) + k] -= s*matrixA[k*order + i].Conjugate(); - } - } - } - } - } - - /// - /// Nonsymmetric reduction to Hessenberg form. - /// - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Order of initial matrix - /// 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. - internal 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; - var mm1O = (m - 1)*order; - for (var i = m; i < order; i++) - { - scale += Math.Abs(matrixH[mm1O + i].Real) + Math.Abs(matrixH[mm1O + i].Imaginary); - } - - if (scale != 0.0) - { - // Compute Householder transformation. - var h = 0.0; - for (var i = order - 1; i >= m; i--) - { - ort[i] = matrixH[mm1O + i]/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[mm1O + m] = 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; - var jO = j*order; - for (var i = order - 1; i >= m; i--) - { - f += ort[i].Conjugate()*matrixH[jO + i]; - } - - f = f/h; - for (var i = m; i < order; i++) - { - matrixH[jO + i] -= 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[j*order + i]; - } - - f = f/h; - for (var j = m; j < order; j++) - { - matrixH[j*order + i] -= f*ort[j].Conjugate(); - } - } - - ort[m] = scale*ort[m]; - matrixH[mm1O + m] *= -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--) - { - var mm1O = (m - 1)*order; - var mm1Om = mm1O + m; - if (matrixH[mm1Om] != Complex.Zero && ort[m] != Complex.Zero) - { - var norm = (matrixH[mm1Om].Real*ort[m].Real) + (matrixH[mm1Om].Imaginary*ort[m].Imaginary); - - for (var i = m + 1; i < order; i++) - { - ort[i] = matrixH[mm1O + i]; - } - - 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++) - { - var im1 = i - 1; - var im1O = im1*order; - var im1Oi = im1O + i; - var iO = i*order; - if (matrixH[im1Oi].Imaginary != 0.0) - { - var y = matrixH[im1Oi]/matrixH[im1Oi].Magnitude; - matrixH[im1Oi] = matrixH[im1Oi].Magnitude; - for (var j = i; j < order; j++) - { - matrixH[j*order + i] *= y.Conjugate(); - } - - for (var j = 0; j <= Math.Min(i + 1, order - 1); j++) - { - matrixH[iO + j] *= y; - } - - for (var j = 0; j < order; j++) - { - dataEv[(i*order) + j] *= y; - } - } - } - } - - /// - /// Nonsymmetric reduction from Hessenberg to real Schur form. - /// - /// Data array of the eigenvectors - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Order of initial matrix - /// 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. - internal static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, Complex[] dataEv, Complex[] matrixH, int order) - { - // Initialize - var n = order - 1; - var eps = Precision.DoublePrecision; - - double norm; - Complex 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 lm1 = l - 1; - var lm1O = lm1*order; - var lO = l*order; - var tst1 = Math.Abs(matrixH[lm1O + lm1].Real) + Math.Abs(matrixH[lm1O + lm1].Imaginary) + Math.Abs(matrixH[lO + l].Real) + Math.Abs(matrixH[lO + l].Imaginary); - if (Math.Abs(matrixH[lm1O + l].Real) < eps*tst1) - { - break; - } - - l--; - } - - var nm1 = n - 1; - var nm1O = nm1*order; - var nO = n*order; - var nOn = nO + n; - // Check for convergence - // One root found - if (l == n) - { - matrixH[nOn] += exshift; - vectorV[n] = matrixH[nOn]; - n--; - iter = 0; - } - else - { - // Form shift - Complex s; - if (iter != 10 && iter != 20) - { - s = matrixH[nOn]; - x = matrixH[nO + nm1]*matrixH[nm1O + n].Real; - - if (x.Real != 0.0 || x.Imaginary != 0.0) - { - y = (matrixH[nm1O + nm1] - 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[nm1O + n].Real) + Math.Abs(matrixH[(n - 2)*order + nm1].Real); - } - - for (var i = 0; i <= n; i++) - { - matrixH[i*order + i] -= s; - } - - exshift += s; - iter++; - - // Reduce to triangle (rows) - for (var i = l + 1; i <= n; i++) - { - var im1 = i - 1; - var im1O = im1*order; - var im1Oim1 = im1O + im1; - s = matrixH[im1O + i].Real; - norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real); - x = matrixH[im1Oim1]/norm; - vectorV[i - 1] = x; - matrixH[im1Oim1] = norm; - matrixH[im1O + i] = new Complex(0.0, s.Real/norm); - - for (var j = i; j < order; j++) - { - var jO = j*order; - y = matrixH[jO + im1]; - z = matrixH[jO + i]; - matrixH[jO + im1] = (x.Conjugate()*y) + (matrixH[im1O + i].Imaginary*z); - matrixH[jO + i] = (x*z) - (matrixH[im1O + i].Imaginary*y); - } - } - - s = matrixH[nOn]; - if (s.Imaginary != 0.0) - { - s /= matrixH[nOn].Magnitude; - matrixH[nOn] = matrixH[nOn].Magnitude; - - for (var j = n + 1; j < order; j++) - { - matrixH[j*order + n] *= s.Conjugate(); - } - } - - // Inverse operation (columns). - for (var j = l + 1; j <= n; j++) - { - x = vectorV[j - 1]; - var jO = j*order; - var jm1 = j - 1; - var jm1O = jm1*order; - var jm1Oj = jm1O + j; - for (var i = 0; i <= j; i++) - { - var jm1Oi = jm1O + i; - z = matrixH[jO + i]; - if (i != j) - { - y = matrixH[jm1Oi]; - matrixH[jm1Oi] = (x*y) + (matrixH[jm1O + j].Imaginary*z); - } - else - { - y = matrixH[jm1Oi].Real; - matrixH[jm1Oi] = new Complex((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[jm1O + j].Imaginary*z.Real), matrixH[jm1Oi].Imaginary); - } - - matrixH[jO + i] = (x.Conjugate()*z) - (matrixH[jm1O + j].Imaginary*y); - } - - for (var i = 0; i < order; i++) - { - y = dataEv[((j - 1)*order) + i]; - z = dataEv[(j*order) + i]; - dataEv[jm1O + i] = (x*y) + (matrixH[jm1Oj].Imaginary*z); - dataEv[jO + i] = (x.Conjugate()*z) - (matrixH[jm1Oj].Imaginary*y); - } - } - - if (s.Imaginary != 0.0) - { - for (var i = 0; i <= n; i++) - { - matrixH[nO + i] *= s; - } - - for (var i = 0; i < order; i++) - { - dataEv[nO + 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[j*order + i].Real) + Math.Abs(matrixH[j*order + i].Imaginary)); - } - } - - if (order == 1) - { - return; - } - - if (norm == 0.0) - { - return; - } - - for (n = order - 1; n > 0; n--) - { - var nO = n*order; - var nOn = nO + n; - x = vectorV[n]; - matrixH[nOn] = 1.0; - - for (var i = n - 1; i >= 0; i--) - { - z = 0.0; - for (var j = i + 1; j <= n; j++) - { - z += matrixH[j*order + i]*matrixH[nO + j]; - } - - y = x - vectorV[i]; - if (y.Real == 0.0 && y.Imaginary == 0.0) - { - y = eps*norm; - } - - matrixH[nO + i] = z/y; - - // Overflow control - var tr = Math.Abs(matrixH[nO + i].Real) + Math.Abs(matrixH[nO + i].Imaginary); - if ((eps*tr)*tr > 1) - { - for (var j = i; j <= n; j++) - { - matrixH[nO + j] = matrixH[nO + j]/tr; - } - } - } - } - - // Back transformation to get eigenvectors of original matrix - for (var j = order - 1; j > 0; j--) - { - var jO = j*order; - for (var i = 0; i < order; i++) - { - z = Complex.Zero; - for (var k = 0; k <= j; k++) - { - z += dataEv[(k*order) + i]*matrixH[jO + k]; - } - - dataEv[jO + i] = z; - } - } - } - /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs index 951c4bad..e54487f7 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -99,724 +99,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { } - /// - /// Reduces a complex Hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations. - /// - /// Source matrix to reduce - /// Output: Arrays for internal storage of real parts of eigenvalues - /// Output: Arrays for internal storage of imaginary parts of eigenvalues - /// Output: Arrays that contains further information about the transformations. - /// Order of initial matrix - /// 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. - internal 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*order + 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[k*order + i].Real) + Math.Abs(matrixA[k*order + i].Imaginary); - } - - if (scale == 0.0f) - { - tau[i - 1] = Complex32.One; - e[i] = 0.0f; - } - else - { - for (var k = 0; k < i; k++) - { - matrixA[k*order + i] /= scale; - h += matrixA[k*order + i].MagnitudeSquared; - } - - Complex32 g = (float) Math.Sqrt(h); - e[i] = scale*g.Real; - - Complex32 temp; - var im1Oi = (i - 1)*order + i; - var f = matrixA[im1Oi]; - if (f.Magnitude != 0.0f) - { - temp = -(matrixA[im1Oi].Conjugate()*tau[i].Conjugate())/f.Magnitude; - h += f.Magnitude*g.Real; - g = 1.0f + (g/f.Magnitude); - matrixA[im1Oi] *= g; - } - else - { - temp = -tau[i].Conjugate(); - matrixA[im1Oi] = g; - } - - if ((f.Magnitude == 0.0f) || (i != 1)) - { - f = Complex32.Zero; - for (var j = 0; j < i; j++) - { - var tmp = Complex32.Zero; - var jO = j*order; - // Form element of A*U. - for (var k = 0; k <= j; k++) - { - tmp += matrixA[k*order + j]*matrixA[k*order + i].Conjugate(); - } - - for (var k = j + 1; k <= i - 1; k++) - { - tmp += matrixA[jO + k].Conjugate()*matrixA[k*order + i].Conjugate(); - } - - // Form element of P - tau[j] = tmp/h; - f += (tmp/h)*matrixA[jO + i]; - } - - hh = f.Real/(h + h); - - // Form the reduced A. - for (var j = 0; j < i; j++) - { - f = matrixA[j*order + i].Conjugate(); - g = tau[j] - (hh*f); - tau[j] = g.Conjugate(); - - for (var k = 0; k <= j; k++) - { - matrixA[k*order + j] -= (f*tau[k]) + (g*matrixA[k*order + i]); - } - } - } - - for (var k = 0; k < i; k++) - { - matrixA[k*order + i] *= scale; - } - - tau[i - 1] = temp.Conjugate(); - } - - hh = d[i]; - d[i] = matrixA[i*order + i].Real; - matrixA[i*order + i] = new Complex32(hh, scale*(float) Math.Sqrt(h)); - } - - hh = d[0]; - d[0] = matrixA[0].Real; - matrixA[0] = hh; - e[0] = 0.0f; - } - - /// - /// Symmetric tridiagonal QL algorithm. - /// - /// Data array of matrix V (eigenvectors) - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - /// - internal 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.DoublePrecision; - 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 NonConvergenceException(); - } - } 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; - } - } - } - } - - /// - /// Determines eigenvectors by undoing the symmetric tridiagonalize transformation - /// - /// Data array of matrix V (eigenvectors) - /// Previously tridiagonalized matrix by . - /// Contains further information about the transformations - /// Input matrix order - /// 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. - internal 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*order + 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[k*order + i]; - } - - s = (s/h)/h; - - for (var k = 0; k < i; k++) - { - dataEv[(j*order) + k] -= s*matrixA[k*order + i].Conjugate(); - } - } - } - } - } - - /// - /// Nonsymmetric reduction to Hessenberg form. - /// - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Order of initial matrix - /// 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. - internal 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; - var mm1O = (m - 1)*order; - for (var i = m; i < order; i++) - { - scale += Math.Abs(matrixH[mm1O + i].Real) + Math.Abs(matrixH[mm1O + i].Imaginary); - } - - if (scale != 0.0f) - { - // Compute Householder transformation. - var h = 0.0f; - for (var i = order - 1; i >= m; i--) - { - ort[i] = matrixH[mm1O + i]/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[mm1O + m] = 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; - var jO = j*order; - for (var i = order - 1; i >= m; i--) - { - f += ort[i].Conjugate()*matrixH[jO + i]; - } - - f = f/h; - for (var i = m; i < order; i++) - { - matrixH[jO + i] -= 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[j*order + i]; - } - - f = f/h; - for (var j = m; j < order; j++) - { - matrixH[j*order + i] -= f*ort[j].Conjugate(); - } - } - - ort[m] = scale*ort[m]; - matrixH[mm1O + m] *= -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--) - { - var mm1O = (m - 1)*order; - var mm1Om = mm1O + m; - if (matrixH[mm1Om] != Complex32.Zero && ort[m] != Complex32.Zero) - { - var norm = (matrixH[mm1Om].Real*ort[m].Real) + (matrixH[mm1Om].Imaginary*ort[m].Imaginary); - - for (var i = m + 1; i < order; i++) - { - ort[i] = matrixH[mm1O + i]; - } - - 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++) - { - var im1 = i - 1; - var im1O = im1*order; - var im1Oi = im1O + i; - var iO = i*order; - if (matrixH[im1Oi].Imaginary != 0.0f) - { - var y = matrixH[im1Oi]/matrixH[im1Oi].Magnitude; - matrixH[im1Oi] = matrixH[im1Oi].Magnitude; - for (var j = i; j < order; j++) - { - matrixH[j*order + i] *= y.Conjugate(); - } - - for (var j = 0; j <= Math.Min(i + 1, order - 1); j++) - { - matrixH[iO + j] *= y; - } - - for (var j = 0; j < order; j++) - { - dataEv[(i*order) + j] *= y; - } - } - } - } - - /// - /// Nonsymmetric reduction from Hessenberg to real Schur form. - /// - /// Data array of the eigenvectors - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Order of initial matrix - /// 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. - internal static void NonsymmetricReduceHessenberToRealSchur(Complex32[] vectorV, Complex32[] dataEv, Complex32[] matrixH, int order) - { - // Initialize - var n = order - 1; - var eps = (float) Precision.SinglePrecision; - - float norm; - Complex32 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 lm1 = l - 1; - var lm1O = lm1*order; - var lO = l*order; - var tst1 = Math.Abs(matrixH[lm1O + lm1].Real) + Math.Abs(matrixH[lm1O + lm1].Imaginary) + Math.Abs(matrixH[lO + l].Real) + Math.Abs(matrixH[lO + l].Imaginary); - if (Math.Abs(matrixH[lm1O + l].Real) < eps*tst1) - { - break; - } - - l--; - } - - var nm1 = n - 1; - var nm1O = nm1*order; - var nO = n*order; - var nOn = nO + n; - // Check for convergence - // One root found - if (l == n) - { - matrixH[nOn] += exshift; - vectorV[n] = matrixH[nOn]; - n--; - iter = 0; - } - else - { - // Form shift - Complex32 s; - if (iter != 10 && iter != 20) - { - s = matrixH[nOn]; - x = matrixH[nO + nm1]*matrixH[nm1O + n].Real; - - if (x.Real != 0.0f || x.Imaginary != 0.0f) - { - y = (matrixH[nm1O + nm1] - s)/2.0f; - z = ((y*y) + x).SquareRoot(); - if ((y.Real*z.Real) + (y.Imaginary*z.Imaginary) < 0.0) - { - z *= -1.0f; - } - - x /= y + z; - s = s - x; - } - } - else - { - // Form exceptional shift - s = Math.Abs(matrixH[nm1O + n].Real) + Math.Abs(matrixH[(n - 2)*order + nm1].Real); - } - - for (var i = 0; i <= n; i++) - { - matrixH[i*order + i] -= s; - } - - exshift += s; - iter++; - - // Reduce to triangle (rows) - for (var i = l + 1; i <= n; i++) - { - var im1 = i - 1; - var im1O = im1*order; - var im1Oim1 = im1O + im1; - s = matrixH[im1O + i].Real; - norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real); - x = matrixH[im1Oim1]/norm; - vectorV[i - 1] = x; - matrixH[im1Oim1] = norm; - matrixH[im1O + i] = new Complex32(0.0f, s.Real/norm); - - for (var j = i; j < order; j++) - { - var jO = j*order; - y = matrixH[jO + im1]; - z = matrixH[jO + i]; - matrixH[jO + im1] = (x.Conjugate()*y) + (matrixH[im1O + i].Imaginary*z); - matrixH[jO + i] = (x*z) - (matrixH[im1O + i].Imaginary*y); - } - } - - s = matrixH[nOn]; - if (s.Imaginary != 0.0f) - { - s /= matrixH[nOn].Magnitude; - matrixH[nOn] = matrixH[nOn].Magnitude; - - for (var j = n + 1; j < order; j++) - { - matrixH[j*order + n] *= s.Conjugate(); - } - } - - // Inverse operation (columns). - for (var j = l + 1; j <= n; j++) - { - x = vectorV[j - 1]; - var jO = j*order; - var jm1 = j - 1; - var jm1O = jm1*order; - var jm1Oj = jm1O + j; - for (var i = 0; i <= j; i++) - { - var jm1Oi = jm1O + i; - z = matrixH[jO + i]; - if (i != j) - { - y = matrixH[jm1Oi]; - matrixH[jm1Oi] = (x*y) + (matrixH[jm1O + j].Imaginary*z); - } - else - { - y = matrixH[jm1Oi].Real; - matrixH[jm1Oi] = new Complex32((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[jm1O + j].Imaginary*z.Real), matrixH[jm1Oi].Imaginary); - } - - matrixH[jO + i] = (x.Conjugate()*z) - (matrixH[jm1O + j].Imaginary*y); - } - - for (var i = 0; i < order; i++) - { - y = dataEv[((j - 1)*order) + i]; - z = dataEv[(j*order) + i]; - dataEv[jm1O + i] = (x*y) + (matrixH[jm1Oj].Imaginary*z); - dataEv[jO + i] = (x.Conjugate()*z) - (matrixH[jm1Oj].Imaginary*y); - } - } - - if (s.Imaginary != 0.0f) - { - for (var i = 0; i <= n; i++) - { - matrixH[nO + i] *= s; - } - - for (var i = 0; i < order; i++) - { - dataEv[nO + 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[j*order + i].Real) + Math.Abs(matrixH[j*order + i].Imaginary)); - } - } - - if (order == 1) - { - return; - } - - if (norm == 0.0) - { - return; - } - - for (n = order - 1; n > 0; n--) - { - var nO = n*order; - var nOn = nO + n; - x = vectorV[n]; - matrixH[nOn] = 1.0f; - - for (var i = n - 1; i >= 0; i--) - { - z = 0.0f; - for (var j = i + 1; j <= n; j++) - { - z += matrixH[j*order + i]*matrixH[nO + j]; - } - - y = x - vectorV[i]; - if (y.Real == 0.0f && y.Imaginary == 0.0f) - { - y = eps*norm; - } - - matrixH[nO + i] = z/y; - - // Overflow control - var tr = Math.Abs(matrixH[nO + i].Real) + Math.Abs(matrixH[nO + i].Imaginary); - if ((eps*tr)*tr > 1) - { - for (var j = i; j <= n; j++) - { - matrixH[nO + j] = matrixH[nO + j]/tr; - } - } - } - } - - // Back transformation to get eigenvectors of original matrix - for (var j = order - 1; j > 0; j--) - { - var jO = j*order; - for (var i = 0; i < order; i++) - { - z = Complex32.Zero; - for (var k = 0; k <= j; k++) - { - z += dataEv[(k*order) + i]*matrixH[jO + k]; - } - - dataEv[jO + i] = z; - } - } - } - /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs index a41a8b02..b61da00c 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -99,1016 +99,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { } - /// - /// Symmetric Householder reduction to tridiagonal form. - /// - /// Data array of matrix V (eigenvectors) - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - internal static void SymmetricTridiagonalize(double[] a, double[] d, double[] e, int order) - { - // 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(d[k]); - } - - if (scale == 0.0) - { - e[i] = d[i - 1]; - for (var j = 0; j < i; j++) - { - d[j] = a[(j*order) + i - 1]; - a[(j*order) + i] = 0.0; - a[(i*order) + j] = 0.0; - } - } - else - { - // Generate Householder vector. - for (var k = 0; k < i; k++) - { - d[k] /= scale; - h += d[k]*d[k]; - } - - var f = d[i - 1]; - var g = Math.Sqrt(h); - if (f > 0) - { - g = -g; - } - - e[i] = scale*g; - h = h - (f*g); - d[i - 1] = f - g; - - for (var j = 0; j < i; j++) - { - e[j] = 0.0; - } - - // Apply similarity transformation to remaining columns. - for (var j = 0; j < i; j++) - { - f = d[j]; - a[(i*order) + j] = f; - g = e[j] + (a[(j*order) + j]*f); - - for (var k = j + 1; k <= i - 1; k++) - { - g += a[(j*order) + k]*d[k]; - e[k] += a[(j*order) + k]*f; - } - - e[j] = g; - } - - f = 0.0; - - for (var j = 0; j < i; j++) - { - e[j] /= h; - f += e[j]*d[j]; - } - - var hh = f/(h + h); - - for (var j = 0; j < i; j++) - { - e[j] -= hh*d[j]; - } - - for (var j = 0; j < i; j++) - { - f = d[j]; - g = e[j]; - - for (var k = j; k <= i - 1; k++) - { - a[(j*order) + k] -= (f*e[k]) + (g*d[k]); - } - - d[j] = a[(j*order) + i - 1]; - a[(j*order) + i] = 0.0; - } - } - - d[i] = h; - } - - // Accumulate transformations. - for (var i = 0; i < order - 1; i++) - { - a[(i*order) + order - 1] = a[(i*order) + i]; - a[(i*order) + i] = 1.0; - var h = d[i + 1]; - if (h != 0.0) - { - for (var k = 0; k <= i; k++) - { - d[k] = a[((i + 1)*order) + k]/h; - } - - for (var j = 0; j <= i; j++) - { - var g = 0.0; - for (var k = 0; k <= i; k++) - { - g += a[((i + 1)*order) + k]*a[(j*order) + k]; - } - - for (var k = 0; k <= i; k++) - { - a[(j*order) + k] -= g*d[k]; - } - } - } - - for (var k = 0; k <= i; k++) - { - a[((i + 1)*order) + k] = 0.0; - } - } - - for (var j = 0; j < order; j++) - { - d[j] = a[(j*order) + order - 1]; - a[(j*order) + order - 1] = 0.0; - } - - a[(order*order) - 1] = 1.0; - e[0] = 0.0; - } - - /// - /// Symmetric tridiagonal QL algorithm. - /// - /// Data array of matrix V (eigenvectors) - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - /// - internal static void SymmetricDiagonalize(double[] a, 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.DoublePrecision; - 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 = a[((i + 1)*order) + k]; - a[((i + 1)*order) + k] = (s*a[(i*order) + k]) + (c*h); - a[(i*order) + k] = (c*a[(i*order) + k]) - (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 NonConvergenceException(); - } - } 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 = a[(i*order) + j]; - a[(i*order) + j] = a[(k*order) + j]; - a[(k*order) + j] = p; - } - } - } - } - - /// - /// Nonsymmetric reduction to Hessenberg form. - /// - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Order of initial matrix - /// 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. - internal static void NonsymmetricReduceToHessenberg(double[] a, double[] matrixH, int order) - { - var ort = new double[order]; - var high = order - 1; - for (var m = 1; m <= high - 1; m++) - { - var mm1 = m - 1; - var mm1O = mm1*order; - // Scale column. - var scale = 0.0; - for (var i = m; i <= high; i++) - { - scale += Math.Abs(matrixH[mm1O + i]); - } - - if (scale != 0.0) - { - // Compute Householder transformation. - var h = 0.0; - for (var i = high; i >= m; i--) - { - ort[i] = matrixH[mm1O + i]/scale; - h += ort[i]*ort[i]; - } - - var g = 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 (var j = m; j < order; j++) - { - var jO = j*order; - var f = 0.0; - for (var i = order - 1; i >= m; i--) - { - f += ort[i]*matrixH[jO + i]; - } - - f = f/h; - - for (var i = m; i <= high; i++) - { - matrixH[jO + i] -= f*ort[i]; - } - } - - for (var i = 0; i <= high; i++) - { - var f = 0.0; - for (var j = high; j >= m; j--) - { - f += ort[j]*matrixH[j*order + i]; - } - f = f/h; - - for (var j = m; j <= high; j++) - { - matrixH[j*order + i] -= f*ort[j]; - } - } - - ort[m] = scale*ort[m]; - matrixH[mm1O + m] = scale*g; - } - } - - // Accumulate transformations (Algol's ortran). - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - a[(j*order) + i] = i == j ? 1.0 : 0.0; - } - } - - for (var m = high - 1; m >= 1; m--) - { - var mm1 = m - 1; - var mm1O = mm1*order; - var mm1Om = mm1O + m; - if (matrixH[mm1Om] != 0.0) - { - for (var i = m + 1; i <= high; i++) - { - ort[i] = matrixH[mm1O + i]; - } - - for (var j = m; j <= high; j++) - { - var g = 0.0; - var jO = j*order; - for (var i = m; i <= high; i++) - { - g += ort[i]*a[jO + i]; - } - - // Double division avoids possible underflow - g = (g/ort[m])/matrixH[mm1Om]; - - for (var i = m; i <= high; i++) - { - a[jO + i] += g*ort[i]; - } - } - } - } - } - - /// - /// Nonsymmetric reduction from Hessenberg to real Schur form. - /// - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - /// - internal static void NonsymmetricReduceHessenberToRealSchur(double[] a, double[] matrixH, double[] d, double[] e, int order) - { - // Initialize - var n = order - 1; - var eps = Math.Pow(2.0, -52.0); - var exshift = 0.0; - double p = 0, q = 0, r = 0, s = 0, z = 0; - double w, x, y; - - // Store roots isolated by balanc and compute matrix norm - var norm = 0.0; - for (var i = 0; i < order; i++) - { - for (var j = Math.Max(i - 1, 0); j < order; j++) - { - norm = norm + Math.Abs(matrixH[j*order + i]); - } - } - - // 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 lm1 = l - 1; - var lm1O = lm1*order; - s = Math.Abs(matrixH[lm1O + lm1]) + Math.Abs(matrixH[l*order + l]); - - if (s == 0.0) - { - s = norm; - } - - if (Math.Abs(matrixH[lm1O + l]) < eps*s) - { - break; - } - - l--; - } - - // Check for convergence - // One root found - if (l == n) - { - var index = n*order + n; - matrixH[index] += exshift; - d[n] = matrixH[index]; - e[n] = 0.0; - n--; - iter = 0; - - // Two roots found - } - else if (l == n - 1) - { - var nO = n*order; - var nm1 = n - 1; - var nm1O = nm1*order; - var nOn = nO + n; - - w = matrixH[nm1O + n]*matrixH[nO + nm1]; - p = (matrixH[nm1O + nm1] - matrixH[nOn])/2.0; - q = (p*p) + w; - z = Math.Sqrt(Math.Abs(q)); - - matrixH[nOn] += exshift; - matrixH[nm1O + nm1] += exshift; - x = matrixH[nOn]; - - // Real pair - if (q >= 0) - { - if (p >= 0) - { - z = p + z; - } - else - { - z = p - z; - } - - d[nm1] = x + z; - - d[n] = d[nm1]; - if (z != 0.0) - { - d[n] = x - (w/z); - } - - e[n - 1] = 0.0; - e[n] = 0.0; - x = matrixH[nm1O + n]; - s = Math.Abs(x) + Math.Abs(z); - p = x/s; - q = z/s; - r = Math.Sqrt((p*p) + (q*q)); - p = p/r; - q = q/r; - - // Row modification - for (var j = n - 1; j < order; j++) - { - var jO = j*order; - var jOn = jO + n; - z = matrixH[jO + nm1]; - matrixH[jO + nm1] = (q*z) + (p*matrixH[jOn]); - matrixH[jOn] = (q*matrixH[jOn]) - (p*z); - } - - // Column modification - for (var i = 0; i <= n; i++) - { - var nOi = nO + i; - z = matrixH[nm1O + i]; - matrixH[nm1O + i] = (q*z) + (p*matrixH[nOi]); - matrixH[nOi] = (q*matrixH[nOi]) - (p*z); - } - - // Accumulate transformations - for (var i = 0; i < order; i++) - { - var nOi = nO + i; - z = a[nm1O + i]; - a[nm1O + i] = (q*z) + (p*a[nOi]); - a[nOi] = (q*a[nOi]) - (p*z); - } - - // Complex pair - } - else - { - d[n - 1] = x + p; - d[n] = x + p; - e[n - 1] = z; - e[n] = -z; - } - - n = n - 2; - iter = 0; - - // No convergence yet - } - else - { - var nO = n*order; - var nm1 = n - 1; - var nm1O = nm1*order; - var nOn = nO + n; - - // Form shift - x = matrixH[nOn]; - y = 0.0; - w = 0.0; - if (l < n) - { - y = matrixH[nm1O + nm1]; - w = matrixH[nm1O + n]*matrixH[nO + nm1]; - } - - // Wilkinson's original ad hoc shift - if (iter == 10) - { - exshift += x; - for (var i = 0; i <= n; i++) - { - matrixH[i*order + i] -= x; - } - - s = Math.Abs(matrixH[nm1O + n]) + Math.Abs(matrixH[(n - 2)*order + nm1]); - x = y = 0.75*s; - w = (-0.4375)*s*s; - } - - // MATLAB's new ad hoc shift - if (iter == 30) - { - s = (y - x)/2.0; - s = (s*s) + w; - if (s > 0) - { - s = Math.Sqrt(s); - if (y < x) - { - s = -s; - } - - s = x - (w/(((y - x)/2.0) + s)); - for (var i = 0; i <= n; i++) - { - matrixH[i*order + i] -= s; - } - - exshift += s; - x = y = w = 0.964; - } - } - - iter = iter + 1; - if (iter >= 30*order) - { - throw new NonConvergenceException(); - } - - // Look for two consecutive small sub-diagonal elements - var m = n - 2; - while (m >= l) - { - var mp1 = m + 1; - var mm1 = m - 1; - var mO = m*order; - var mp1O = mp1*order; - var mm1O = mm1*order; - - z = matrixH[mO + m]; - r = x - z; - s = y - z; - p = (((r*s) - w)/matrixH[mO + mp1]) + matrixH[mp1O + m]; - q = matrixH[mp1O + mp1] - z - r - s; - r = matrixH[mp1O + (m + 2)]; - s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); - p = p/s; - q = q/s; - r = r/s; - - if (m == l) - { - break; - } - - if (Math.Abs(matrixH[mm1O + m])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[mm1O + mm1]) + Math.Abs(z) + Math.Abs(matrixH[mp1O + mp1])))) - { - break; - } - - m--; - } - - var mp2 = m + 2; - for (var i = mp2; i <= n; i++) - { - matrixH[(i - 2)*order + i] = 0.0; - if (i > mp2) - { - matrixH[(i - 3)*order + i] = 0.0; - } - } - - // Double QR step involving rows l:n and columns m:n - for (var k = m; k <= n - 1; k++) - { - var notlast = k != n - 1; - var kO = k*order; - var km1 = k - 1; - var kp1 = k + 1; - var kp2 = k + 2; - var kp1O = kp1*order; - var kp2O = kp2*order; - var km1O = km1*order; - if (k != m) - { - p = matrixH[km1O + k]; - q = matrixH[km1O + kp1]; - r = notlast ? matrixH[km1O + kp2] : 0.0; - x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); - if (x == 0.0) - { - continue; - } - - p = p/x; - q = q/x; - r = r/x; - } - - s = Math.Sqrt((p*p) + (q*q) + (r*r)); - - if (p < 0) - { - s = -s; - } - - if (s != 0.0) - { - if (k != m) - { - matrixH[km1O + k] = (-s)*x; - } - else if (l != m) - { - matrixH[km1O + k] = -matrixH[km1O + k]; - } - - p = p + s; - x = p/s; - y = q/s; - z = r/s; - q = q/p; - r = r/p; - - // Row modification - for (var j = k; j < order; j++) - { - var jO = j*order; - var jOk = jO + k; - var jOkp1 = jO + kp1; - var jOkp2 = jO + kp2; - p = matrixH[jOk] + (q*matrixH[jOkp1]); - if (notlast) - { - p = p + (r*matrixH[jOkp2]); - matrixH[jOkp2] -= (p*z); - } - - matrixH[jOk] -= (p*x); - matrixH[jOkp1] -= (p*y); - } - - // Column modification - for (var i = 0; i <= Math.Min(n, k + 3); i++) - { - p = (x*matrixH[kO + i]) + (y*matrixH[kp1O + i]); - - if (notlast) - { - p = p + (z*matrixH[kp2O + i]); - matrixH[kp2O + i] -= (p*r); - } - - matrixH[kO + i] -= p; - matrixH[kp1O + i] -= (p*q); - } - - // Accumulate transformations - for (var i = 0; i < order; i++) - { - p = (x*a[kO + i]) + (y*a[kp1O + i]); - - if (notlast) - { - p = p + (z*a[kp2O + i]); - a[kp2O + i] -= p*r; - } - - a[kO + i] -= p; - a[kp1O + i] -= p*q; - } - } // (s != 0) - } // k loop - } // check convergence - } // while (n >= low) - - // Backsubstitute to find vectors of upper triangular form - if (norm == 0.0) - { - return; - } - - for (n = order - 1; n >= 0; n--) - { - var nO = n*order; - var nm1 = n - 1; - var nm1O = nm1*order; - - p = d[n]; - q = e[n]; - - - // Real vector - double t; - if (q == 0.0) - { - var l = n; - matrixH[nO + n] = 1.0; - for (var i = n - 1; i >= 0; i--) - { - var ip1 = i + 1; - var iO = i*order; - var ip1O = ip1*order; - - w = matrixH[iO + i] - p; - r = 0.0; - for (var j = l; j <= n; j++) - { - r = r + (matrixH[j*order + i]*matrixH[nO + j]); - } - - if (e[i] < 0.0) - { - z = w; - s = r; - } - else - { - l = i; - if (e[i] == 0.0) - { - if (w != 0.0) - { - matrixH[nO + i] = (-r)/w; - } - else - { - matrixH[nO + i] = (-r)/(eps*norm); - } - - // Solve real equations - } - else - { - x = matrixH[ip1O + i]; - y = matrixH[iO + ip1]; - q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]); - t = ((x*s) - (z*r))/q; - matrixH[nO + i] = t; - if (Math.Abs(x) > Math.Abs(z)) - { - matrixH[nO + ip1] = (-r - (w*t))/x; - } - else - { - matrixH[nO + ip1] = (-s - (y*t))/z; - } - } - - // Overflow control - t = Math.Abs(matrixH[nO + i]); - if ((eps*t)*t > 1) - { - for (var j = i; j <= n; j++) - { - matrixH[nO + j] /= t; - } - } - } - } - - // Complex vector - } - else if (q < 0) - { - var l = n - 1; - - // Last vector component imaginary so matrix is triangular - if (Math.Abs(matrixH[nm1O + n]) > Math.Abs(matrixH[nO + nm1])) - { - matrixH[nm1O + nm1] = q/matrixH[nm1O + n]; - matrixH[nO + nm1] = (-(matrixH[nO + n] - p))/matrixH[nm1O + n]; - } - else - { - var res = Cdiv(0.0, -matrixH[nO + nm1], matrixH[nm1O + nm1] - p, q); - matrixH[nm1O + nm1] = res.Real; - matrixH[nO + nm1] = res.Imaginary; - } - - matrixH[nm1O + n] = 0.0; - matrixH[nO + n] = 1.0; - for (var i = n - 2; i >= 0; i--) - { - var ip1 = i + 1; - var iO = i*order; - var ip1O = ip1*order; - var ra = 0.0; - var sa = 0.0; - for (var j = l; j <= n; j++) - { - var jO = j*order; - var jOi = jO + i; - ra = ra + (matrixH[jOi]*matrixH[nm1O + j]); - sa = sa + (matrixH[jOi]*matrixH[nO + j]); - } - - w = matrixH[iO + i] - p; - - if (e[i] < 0.0) - { - z = w; - r = ra; - s = sa; - } - else - { - l = i; - if (e[i] == 0.0) - { - var res = Cdiv(-ra, -sa, w, q); - matrixH[nm1O + i] = res.Real; - matrixH[nO + i] = res.Imaginary; - } - else - { - // Solve complex equations - x = matrixH[ip1O + i]; - y = matrixH[iO + ip1]; - - var vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q); - var vi = (d[i] - p)*2.0*q; - if ((vr == 0.0) && (vi == 0.0)) - { - vr = eps*norm*(Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); - } - - var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi); - matrixH[nm1O + i] = res.Real; - matrixH[nO + i] = res.Imaginary; - if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q))) - { - matrixH[nm1O + ip1] = (-ra - (w*matrixH[nm1O + i]) + (q*matrixH[nO + i]))/x; - matrixH[nO + ip1] = (-sa - (w*matrixH[nO + i]) - (q*matrixH[nm1O + i]))/x; - } - else - { - res = Cdiv(-r - (y*matrixH[nm1O + i]), -s - (y*matrixH[nO + i]), z, q); - matrixH[nm1O + ip1] = res.Real; - matrixH[nO + ip1] = res.Imaginary; - } - } - - // Overflow control - t = Math.Max(Math.Abs(matrixH[nm1O + i]), Math.Abs(matrixH[nO + i])); - if ((eps*t)*t > 1) - { - for (var j = i; j <= n; j++) - { - matrixH[nm1O + j] /= t; - matrixH[nO + j] /= t; - } - } - } - } - } - } - - // Back transformation to get eigenvectors of original matrix - for (var j = order - 1; j >= 0; j--) - { - var jO = j*order; - for (var i = 0; i < order; i++) - { - z = 0.0; - for (var k = 0; k <= j; k++) - { - z = z + (a[k*order + i]*matrixH[jO + k]); - } - - a[jO + i] = z; - } - } - } - - /// - /// Complex scalar division X/Y. - /// - /// Real part of X - /// Imaginary part of X - /// Real part of Y - /// Imaginary part of Y - /// Division result as a number. - static Complex Cdiv(double xreal, double ximag, double yreal, double yimag) - { - if (Math.Abs(yimag) < Math.Abs(yreal)) - { - return new Complex((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal)))); - } - - return new Complex((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag)))); - } - /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs index f684f2ab..2dda05a9 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -99,1016 +99,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { } - /// - /// Symmetric Householder reduction to tridiagonal form. - /// - /// Data array of matrix V (eigenvectors) - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - internal static void SymmetricTridiagonalize(float[] a, float[] d, float[] e, int order) - { - // 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(d[k]); - } - - if (scale == 0.0f) - { - e[i] = d[i - 1]; - for (var j = 0; j < i; j++) - { - d[j] = a[(j*order) + i - 1]; - a[(j*order) + i] = 0.0f; - a[(i*order) + j] = 0.0f; - } - } - else - { - // Generate Householder vector. - for (var k = 0; k < i; k++) - { - d[k] /= scale; - h += d[k]*d[k]; - } - - var f = d[i - 1]; - var g = (float) Math.Sqrt(h); - if (f > 0) - { - g = -g; - } - - e[i] = scale*g; - h = h - (f*g); - d[i - 1] = f - g; - - for (var j = 0; j < i; j++) - { - e[j] = 0.0f; - } - - // Apply similarity transformation to remaining columns. - for (var j = 0; j < i; j++) - { - f = d[j]; - a[(i*order) + j] = f; - g = e[j] + (a[(j*order) + j]*f); - - for (var k = j + 1; k <= i - 1; k++) - { - g += a[(j*order) + k]*d[k]; - e[k] += a[(j*order) + k]*f; - } - - e[j] = g; - } - - f = 0.0f; - - for (var j = 0; j < i; j++) - { - e[j] /= h; - f += e[j]*d[j]; - } - - var hh = f/(h + h); - - for (var j = 0; j < i; j++) - { - e[j] -= hh*d[j]; - } - - for (var j = 0; j < i; j++) - { - f = d[j]; - g = e[j]; - - for (var k = j; k <= i - 1; k++) - { - a[(j*order) + k] -= (f*e[k]) + (g*d[k]); - } - - d[j] = a[(j*order) + i - 1]; - a[(j*order) + i] = 0.0f; - } - } - - d[i] = h; - } - - // Accumulate transformations. - for (var i = 0; i < order - 1; i++) - { - a[(i*order) + order - 1] = a[(i*order) + i]; - a[(i*order) + i] = 1.0f; - var h = d[i + 1]; - if (h != 0.0f) - { - for (var k = 0; k <= i; k++) - { - d[k] = a[((i + 1)*order) + k]/h; - } - - for (var j = 0; j <= i; j++) - { - var g = 0.0f; - for (var k = 0; k <= i; k++) - { - g += a[((i + 1)*order) + k]*a[(j*order) + k]; - } - - for (var k = 0; k <= i; k++) - { - a[(j*order) + k] -= g*d[k]; - } - } - } - - for (var k = 0; k <= i; k++) - { - a[((i + 1)*order) + k] = 0.0f; - } - } - - for (var j = 0; j < order; j++) - { - d[j] = a[(j*order) + order - 1]; - a[(j*order) + order - 1] = 0.0f; - } - - a[(order*order) - 1] = 1.0f; - e[0] = 0.0f; - } - - /// - /// Symmetric tridiagonal QL algorithm. - /// - /// Data array of matrix V (eigenvectors) - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - /// - internal static void SymmetricDiagonalize(float[] a, 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.SinglePrecision; - 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 = a[((i + 1)*order) + k]; - a[((i + 1)*order) + k] = (s*a[(i*order) + k]) + (c*h); - a[(i*order) + k] = (c*a[(i*order) + k]) - (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 NonConvergenceException(); - } - } 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 = a[(i*order) + j]; - a[(i*order) + j] = a[(k*order) + j]; - a[(k*order) + j] = p; - } - } - } - } - - /// - /// Nonsymmetric reduction to Hessenberg form. - /// - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Order of initial matrix - /// 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. - internal static void NonsymmetricReduceToHessenberg(float[] a, float[] matrixH, int order) - { - var ort = new float[order]; - var high = order - 1; - for (var m = 1; m <= high - 1; m++) - { - var mm1 = m - 1; - var mm1O = mm1*order; - // Scale column. - var scale = 0.0f; - for (var i = m; i <= high; i++) - { - scale += Math.Abs(matrixH[mm1O + i]); - } - - if (scale != 0.0f) - { - // Compute Householder transformation. - var h = 0.0f; - for (var i = high; i >= m; i--) - { - ort[i] = matrixH[mm1O + i]/scale; - h += ort[i]*ort[i]; - } - - var g = (float) 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 (var j = m; j < order; j++) - { - var jO = j*order; - var f = 0.0f; - for (var i = order - 1; i >= m; i--) - { - f += ort[i]*matrixH[jO + i]; - } - - f = f/h; - - for (var i = m; i <= high; i++) - { - matrixH[jO + i] -= f*ort[i]; - } - } - - for (var i = 0; i <= high; i++) - { - var f = 0.0f; - for (var j = high; j >= m; j--) - { - f += ort[j]*matrixH[j*order + i]; - } - f = f/h; - - for (var j = m; j <= high; j++) - { - matrixH[j*order + i] -= f*ort[j]; - } - } - - ort[m] = scale*ort[m]; - matrixH[mm1O + m] = scale*g; - } - } - - // Accumulate transformations (Algol's ortran). - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - a[(j*order) + i] = i == j ? 1.0f : 0.0f; - } - } - - for (var m = high - 1; m >= 1; m--) - { - var mm1 = m - 1; - var mm1O = mm1*order; - var mm1Om = mm1O + m; - if (matrixH[mm1Om] != 0.0) - { - for (var i = m + 1; i <= high; i++) - { - ort[i] = matrixH[mm1O + i]; - } - - for (var j = m; j <= high; j++) - { - var g = 0.0f; - var jO = j*order; - for (var i = m; i <= high; i++) - { - g += ort[i]*a[jO + i]; - } - - // Double division avoids possible underflow - g = (g/ort[m])/matrixH[mm1Om]; - - for (var i = m; i <= high; i++) - { - a[jO + i] += g*ort[i]; - } - } - } - } - } - - /// - /// Nonsymmetric reduction from Hessenberg to real Schur form. - /// - /// Data array of matrix V (eigenvectors) - /// Array for internal storage of nonsymmetric Hessenberg form. - /// Arrays for internal storage of real parts of eigenvalues - /// Arrays for internal storage of imaginary parts of eigenvalues - /// Order of initial matrix - /// 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. - /// - internal static void NonsymmetricReduceHessenberToRealSchur(float[] a, float[] matrixH, float[] d, float[] e, int order) - { - // Initialize - var n = order - 1; - var eps = (float) Precision.SinglePrecision; - var exshift = 0.0f; - float p = 0, q = 0, r = 0, s = 0, z = 0; - float w, x, y; - - // Store roots isolated by balanc and compute matrix norm - var norm = 0.0f; - for (var i = 0; i < order; i++) - { - for (var j = Math.Max(i - 1, 0); j < order; j++) - { - norm = norm + Math.Abs(matrixH[j*order + i]); - } - } - - // 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 lm1 = l - 1; - var lm1O = lm1*order; - s = Math.Abs(matrixH[lm1O + lm1]) + Math.Abs(matrixH[l*order + l]); - - if (s == 0.0) - { - s = norm; - } - - if (Math.Abs(matrixH[lm1O + l]) < eps*s) - { - break; - } - - l--; - } - - // Check for convergence - // One root found - if (l == n) - { - var index = n*order + n; - matrixH[index] += exshift; - d[n] = matrixH[index]; - e[n] = 0.0f; - n--; - iter = 0; - - // Two roots found - } - else if (l == n - 1) - { - var nO = n*order; - var nm1 = n - 1; - var nm1O = nm1*order; - var nOn = nO + n; - - w = matrixH[nm1O + n]*matrixH[nO + nm1]; - p = (matrixH[nm1O + nm1] - matrixH[nOn])/2.0f; - q = (p*p) + w; - z = (float) Math.Sqrt(Math.Abs(q)); - - matrixH[nOn] += exshift; - matrixH[nm1O + nm1] += exshift; - x = matrixH[nOn]; - - // Real pair - if (q >= 0) - { - if (p >= 0) - { - z = p + z; - } - else - { - z = p - z; - } - - d[nm1] = x + z; - - d[n] = d[nm1]; - if (z != 0.0) - { - d[n] = x - (w/z); - } - - e[n - 1] = 0.0f; - e[n] = 0.0f; - x = matrixH[nm1O + n]; - s = Math.Abs(x) + Math.Abs(z); - p = x/s; - q = z/s; - r = (float) Math.Sqrt((p*p) + (q*q)); - p = p/r; - q = q/r; - - // Row modification - for (var j = n - 1; j < order; j++) - { - var jO = j*order; - var jOn = jO + n; - z = matrixH[jO + nm1]; - matrixH[jO + nm1] = (q*z) + (p*matrixH[jOn]); - matrixH[jOn] = (q*matrixH[jOn]) - (p*z); - } - - // Column modification - for (var i = 0; i <= n; i++) - { - var nOi = nO + i; - z = matrixH[nm1O + i]; - matrixH[nm1O + i] = (q*z) + (p*matrixH[nOi]); - matrixH[nOi] = (q*matrixH[nOi]) - (p*z); - } - - // Accumulate transformations - for (var i = 0; i < order; i++) - { - var nOi = nO + i; - z = a[nm1O + i]; - a[nm1O + i] = (q*z) + (p*a[nOi]); - a[nOi] = (q*a[nOi]) - (p*z); - } - - // Complex pair - } - else - { - d[n - 1] = x + p; - d[n] = x + p; - e[n - 1] = z; - e[n] = -z; - } - - n = n - 2; - iter = 0; - - // No convergence yet - } - else - { - var nO = n*order; - var nm1 = n - 1; - var nm1O = nm1*order; - var nOn = nO + n; - - // Form shift - x = matrixH[nOn]; - y = 0.0f; - w = 0.0f; - if (l < n) - { - y = matrixH[nm1O + nm1]; - w = matrixH[nm1O + n]*matrixH[nO + nm1]; - } - - // Wilkinson's original ad hoc shift - if (iter == 10) - { - exshift += x; - for (var i = 0; i <= n; i++) - { - matrixH[i*order + i] -= x; - } - - s = Math.Abs(matrixH[nm1O + n]) + Math.Abs(matrixH[(n - 2)*order + nm1]); - x = y = 0.75f*s; - w = (-0.4375f)*s*s; - } - - // MATLAB's new ad hoc shift - if (iter == 30) - { - s = (y - x)/2.0f; - s = (s*s) + w; - if (s > 0) - { - s = (float) Math.Sqrt(s); - if (y < x) - { - s = -s; - } - - s = x - (w/(((y - x)/2.0f) + s)); - for (var i = 0; i <= n; i++) - { - matrixH[i*order + i] -= s; - } - - exshift += s; - x = y = w = 0.964f; - } - } - - iter = iter + 1; - if (iter >= 30*order) - { - throw new NonConvergenceException(); - } - - // Look for two consecutive small sub-diagonal elements - var m = n - 2; - while (m >= l) - { - var mp1 = m + 1; - var mm1 = m - 1; - var mO = m*order; - var mp1O = mp1*order; - var mm1O = mm1*order; - - z = matrixH[mO + m]; - r = x - z; - s = y - z; - p = (((r*s) - w)/matrixH[mO + mp1]) + matrixH[mp1O + m]; - q = matrixH[mp1O + mp1] - z - r - s; - r = matrixH[mp1O + (m + 2)]; - s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); - p = p/s; - q = q/s; - r = r/s; - - if (m == l) - { - break; - } - - if (Math.Abs(matrixH[mm1O + m])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[mm1O + mm1]) + Math.Abs(z) + Math.Abs(matrixH[mp1O + mp1])))) - { - break; - } - - m--; - } - - var mp2 = m + 2; - for (var i = mp2; i <= n; i++) - { - matrixH[(i - 2)*order + i] = 0.0f; - if (i > mp2) - { - matrixH[(i - 3)*order + i] = 0.0f; - } - } - - // Double QR step involving rows l:n and columns m:n - for (var k = m; k <= n - 1; k++) - { - var notlast = k != n - 1; - var kO = k*order; - var km1 = k - 1; - var kp1 = k + 1; - var kp2 = k + 2; - var kp1O = kp1*order; - var kp2O = kp2*order; - var km1O = km1*order; - if (k != m) - { - p = matrixH[km1O + k]; - q = matrixH[km1O + kp1]; - r = notlast ? matrixH[km1O + kp2] : 0.0f; - x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); - if (x == 0.0f) - { - continue; - } - - p = p/x; - q = q/x; - r = r/x; - } - - s = (float) Math.Sqrt((p*p) + (q*q) + (r*r)); - - if (p < 0) - { - s = -s; - } - - if (s != 0.0f) - { - if (k != m) - { - matrixH[km1O + k] = (-s)*x; - } - else if (l != m) - { - matrixH[km1O + k] = -matrixH[km1O + k]; - } - - p = p + s; - x = p/s; - y = q/s; - z = r/s; - q = q/p; - r = r/p; - - // Row modification - for (var j = k; j < order; j++) - { - var jO = j*order; - var jOk = jO + k; - var jOkp1 = jO + kp1; - var jOkp2 = jO + kp2; - p = matrixH[jOk] + (q*matrixH[jOkp1]); - if (notlast) - { - p = p + (r*matrixH[jOkp2]); - matrixH[jOkp2] -= (p*z); - } - - matrixH[jOk] -= (p*x); - matrixH[jOkp1] -= (p*y); - } - - // Column modification - for (var i = 0; i <= Math.Min(n, k + 3); i++) - { - p = (x*matrixH[kO + i]) + (y*matrixH[kp1O + i]); - - if (notlast) - { - p = p + (z*matrixH[kp2O + i]); - matrixH[kp2O + i] -= (p*r); - } - - matrixH[kO + i] -= p; - matrixH[kp1O + i] -= (p*q); - } - - // Accumulate transformations - for (var i = 0; i < order; i++) - { - p = (x*a[kO + i]) + (y*a[kp1O + i]); - - if (notlast) - { - p = p + (z*a[kp2O + i]); - a[kp2O + i] -= p*r; - } - - a[kO + i] -= p; - a[kp1O + i] -= p*q; - } - } // (s != 0) - } // k loop - } // check convergence - } // while (n >= low) - - // Backsubstitute to find vectors of upper triangular form - if (norm == 0.0f) - { - return; - } - - for (n = order - 1; n >= 0; n--) - { - var nO = n*order; - var nm1 = n - 1; - var nm1O = nm1*order; - - p = d[n]; - q = e[n]; - - - // Real vector - float t; - if (q == 0.0f) - { - var l = n; - matrixH[nO + n] = 1.0f; - for (var i = n - 1; i >= 0; i--) - { - var ip1 = i + 1; - var iO = i*order; - var ip1O = ip1*order; - - w = matrixH[iO + i] - p; - r = 0.0f; - for (var j = l; j <= n; j++) - { - r = r + (matrixH[j*order + i]*matrixH[nO + j]); - } - - if (e[i] < 0.0) - { - z = w; - s = r; - } - else - { - l = i; - if (e[i] == 0.0f) - { - if (w != 0.0f) - { - matrixH[nO + i] = (-r)/w; - } - else - { - matrixH[nO + i] = (-r)/(eps*norm); - } - - // Solve real equations - } - else - { - x = matrixH[ip1O + i]; - y = matrixH[iO + ip1]; - q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]); - t = ((x*s) - (z*r))/q; - matrixH[nO + i] = t; - if (Math.Abs(x) > Math.Abs(z)) - { - matrixH[nO + ip1] = (-r - (w*t))/x; - } - else - { - matrixH[nO + ip1] = (-s - (y*t))/z; - } - } - - // Overflow control - t = Math.Abs(matrixH[nO + i]); - if ((eps*t)*t > 1) - { - for (var j = i; j <= n; j++) - { - matrixH[nO + j] /= t; - } - } - } - } - - // Complex vector - } - else if (q < 0) - { - var l = n - 1; - - // Last vector component imaginary so matrix is triangular - if (Math.Abs(matrixH[nm1O + n]) > Math.Abs(matrixH[nO + nm1])) - { - matrixH[nm1O + nm1] = q/matrixH[nm1O + n]; - matrixH[nO + nm1] = (-(matrixH[nO + n] - p))/matrixH[nm1O + n]; - } - else - { - var res = Cdiv(0.0f, -matrixH[nO + nm1], matrixH[nm1O + nm1] - p, q); - matrixH[nm1O + nm1] = res.Real; - matrixH[nO + nm1] = res.Imaginary; - } - - matrixH[nm1O + n] = 0.0f; - matrixH[nO + n] = 1.0f; - for (var i = n - 2; i >= 0; i--) - { - var ip1 = i + 1; - var iO = i*order; - var ip1O = ip1*order; - var ra = 0.0f; - var sa = 0.0f; - for (var j = l; j <= n; j++) - { - var jO = j*order; - var jOi = jO + i; - ra = ra + (matrixH[jOi]*matrixH[nm1O + j]); - sa = sa + (matrixH[jOi]*matrixH[nO + j]); - } - - w = matrixH[iO + i] - p; - - if (e[i] < 0.0) - { - z = w; - r = ra; - s = sa; - } - else - { - l = i; - if (e[i] == 0.0) - { - var res = Cdiv(-ra, -sa, w, q); - matrixH[nm1O + i] = res.Real; - matrixH[nO + i] = res.Imaginary; - } - else - { - // Solve complex equations - x = matrixH[ip1O + i]; - y = matrixH[iO + ip1]; - - var vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q); - var vi = (d[i] - p)*2.0f*q; - if ((vr == 0.0f) && (vi == 0.0f)) - { - vr = eps*norm*(Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); - } - - var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi); - matrixH[nm1O + i] = res.Real; - matrixH[nO + i] = res.Imaginary; - if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q))) - { - matrixH[nm1O + ip1] = (-ra - (w*matrixH[nm1O + i]) + (q*matrixH[nO + i]))/x; - matrixH[nO + ip1] = (-sa - (w*matrixH[nO + i]) - (q*matrixH[nm1O + i]))/x; - } - else - { - res = Cdiv(-r - (y*matrixH[nm1O + i]), -s - (y*matrixH[nO + i]), z, q); - matrixH[nm1O + ip1] = res.Real; - matrixH[nO + ip1] = res.Imaginary; - } - } - - // Overflow control - t = Math.Max(Math.Abs(matrixH[nm1O + i]), Math.Abs(matrixH[nO + i])); - if ((eps*t)*t > 1) - { - for (var j = i; j <= n; j++) - { - matrixH[nm1O + j] /= t; - matrixH[nO + j] /= t; - } - } - } - } - } - } - - // Back transformation to get eigenvectors of original matrix - for (var j = order - 1; j >= 0; j--) - { - var jO = j*order; - for (var i = 0; i < order; i++) - { - z = 0.0f; - for (var k = 0; k <= j; k++) - { - z = z + (a[k*order + i]*matrixH[jO + k]); - } - - a[jO + i] = z; - } - } - } - - /// - /// Complex scalar division X/Y. - /// - /// Real part of X - /// Imaginary part of X - /// Real part of Y - /// Imaginary part of Y - /// Division result as a number. - static Numerics.Complex32 Cdiv(float xreal, float ximag, float yreal, float yimag) - { - if (Math.Abs(yimag) < Math.Abs(yreal)) - { - return new Numerics.Complex32((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal)))); - } - - return new Numerics.Complex32((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag)))); - } - /// /// Solves a system of linear equations, AX = B, with A SVD factorized. /// diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs index 6a7d05f4..17d3a3c7 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,11 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Complex.Factorization; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.Managed { @@ -2488,9 +2487,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed var d = new double[order]; var e = new double[order]; - DenseEvd.SymmetricTridiagonalize(matrixCopy, d, e, tau, order); - DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); - DenseEvd.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); + SymmetricTridiagonalize(matrixCopy, d, e, tau, order); + SymmetricDiagonalize(matrixEv, d, e, order); + SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); for (var i = 0; i < order; i++) { @@ -2499,8 +2498,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed } else { - DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); - DenseEvd.NonsymmetricReduceHessenberToRealSchur(vectorEv, matrixEv, matrixCopy, order); + NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); + NonsymmetricReduceHessenberToRealSchur(vectorEv, matrixEv, matrixCopy, order); } for (var i = 0; i < order; i++) @@ -2509,6 +2508,724 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed } } + /// + /// Reduces a complex Hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations. + /// + /// Source matrix to reduce + /// Output: Arrays for internal storage of real parts of eigenvalues + /// Output: Arrays for internal storage of imaginary parts of eigenvalues + /// Output: Arrays that contains further information about the transformations. + /// Order of initial matrix + /// 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. + internal 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*order + 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[k*order + i].Real) + Math.Abs(matrixA[k*order + i].Imaginary); + } + + if (scale == 0.0) + { + tau[i - 1] = Complex.One; + e[i] = 0.0; + } + else + { + for (var k = 0; k < i; k++) + { + matrixA[k*order + i] /= scale; + h += matrixA[k*order + i].MagnitudeSquared(); + } + + Complex g = Math.Sqrt(h); + e[i] = scale*g.Real; + + Complex temp; + var im1Oi = (i - 1)*order + i; + var f = matrixA[im1Oi]; + if (f.Magnitude != 0) + { + temp = -(matrixA[im1Oi].Conjugate()*tau[i].Conjugate())/f.Magnitude; + h += f.Magnitude*g.Real; + g = 1.0 + (g/f.Magnitude); + matrixA[im1Oi] *= g; + } + else + { + temp = -tau[i].Conjugate(); + matrixA[im1Oi] = g; + } + + if ((f.Magnitude == 0) || (i != 1)) + { + f = Complex.Zero; + for (var j = 0; j < i; j++) + { + var tmp = Complex.Zero; + var jO = j*order; + // Form element of A*U. + for (var k = 0; k <= j; k++) + { + tmp += matrixA[k*order + j]*matrixA[k*order + i].Conjugate(); + } + + for (var k = j + 1; k <= i - 1; k++) + { + tmp += matrixA[jO + k].Conjugate()*matrixA[k*order + i].Conjugate(); + } + + // Form element of P + tau[j] = tmp/h; + f += (tmp/h)*matrixA[jO + i]; + } + + hh = f.Real/(h + h); + + // Form the reduced A. + for (var j = 0; j < i; j++) + { + f = matrixA[j*order + i].Conjugate(); + g = tau[j] - (hh*f); + tau[j] = g.Conjugate(); + + for (var k = 0; k <= j; k++) + { + matrixA[k*order + j] -= (f*tau[k]) + (g*matrixA[k*order + i]); + } + } + } + + for (var k = 0; k < i; k++) + { + matrixA[k*order + i] *= scale; + } + + tau[i - 1] = temp.Conjugate(); + } + + hh = d[i]; + d[i] = matrixA[i*order + i].Real; + matrixA[i*order + i] = new Complex(hh, scale*Math.Sqrt(h)); + } + + hh = d[0]; + d[0] = matrixA[0].Real; + matrixA[0] = hh; + e[0] = 0.0; + } + + /// + /// Symmetric tridiagonal QL algorithm. + /// + /// Data array of matrix V (eigenvectors) + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + /// + internal 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.DoublePrecision; + 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 NonConvergenceException(); + } + } 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; + } + } + } + } + + /// + /// Determines eigenvectors by undoing the symmetric tridiagonalize transformation + /// + /// Data array of matrix V (eigenvectors) + /// Previously tridiagonalized matrix by SymmetricTridiagonalize. + /// Contains further information about the transformations + /// Input matrix order + /// 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. + internal 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*order + 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[k*order + i]; + } + + s = (s/h)/h; + + for (var k = 0; k < i; k++) + { + dataEv[(j*order) + k] -= s*matrixA[k*order + i].Conjugate(); + } + } + } + } + } + + /// + /// Nonsymmetric reduction to Hessenberg form. + /// + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Order of initial matrix + /// 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. + internal 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; + var mm1O = (m - 1)*order; + for (var i = m; i < order; i++) + { + scale += Math.Abs(matrixH[mm1O + i].Real) + Math.Abs(matrixH[mm1O + i].Imaginary); + } + + if (scale != 0.0) + { + // Compute Householder transformation. + var h = 0.0; + for (var i = order - 1; i >= m; i--) + { + ort[i] = matrixH[mm1O + i]/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[mm1O + m] = 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; + var jO = j*order; + for (var i = order - 1; i >= m; i--) + { + f += ort[i].Conjugate()*matrixH[jO + i]; + } + + f = f/h; + for (var i = m; i < order; i++) + { + matrixH[jO + i] -= 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[j*order + i]; + } + + f = f/h; + for (var j = m; j < order; j++) + { + matrixH[j*order + i] -= f*ort[j].Conjugate(); + } + } + + ort[m] = scale*ort[m]; + matrixH[mm1O + m] *= -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--) + { + var mm1O = (m - 1)*order; + var mm1Om = mm1O + m; + if (matrixH[mm1Om] != Complex.Zero && ort[m] != Complex.Zero) + { + var norm = (matrixH[mm1Om].Real*ort[m].Real) + (matrixH[mm1Om].Imaginary*ort[m].Imaginary); + + for (var i = m + 1; i < order; i++) + { + ort[i] = matrixH[mm1O + i]; + } + + 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++) + { + var im1 = i - 1; + var im1O = im1*order; + var im1Oi = im1O + i; + var iO = i*order; + if (matrixH[im1Oi].Imaginary != 0.0) + { + var y = matrixH[im1Oi]/matrixH[im1Oi].Magnitude; + matrixH[im1Oi] = matrixH[im1Oi].Magnitude; + for (var j = i; j < order; j++) + { + matrixH[j*order + i] *= y.Conjugate(); + } + + for (var j = 0; j <= Math.Min(i + 1, order - 1); j++) + { + matrixH[iO + j] *= y; + } + + for (var j = 0; j < order; j++) + { + dataEv[(i*order) + j] *= y; + } + } + } + } + + /// + /// Nonsymmetric reduction from Hessenberg to real Schur form. + /// + /// Data array of the eigenvectors + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Order of initial matrix + /// 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. + internal static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, Complex[] dataEv, Complex[] matrixH, int order) + { + // Initialize + var n = order - 1; + var eps = Precision.DoublePrecision; + + double norm; + Complex 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 lm1 = l - 1; + var lm1O = lm1*order; + var lO = l*order; + var tst1 = Math.Abs(matrixH[lm1O + lm1].Real) + Math.Abs(matrixH[lm1O + lm1].Imaginary) + Math.Abs(matrixH[lO + l].Real) + Math.Abs(matrixH[lO + l].Imaginary); + if (Math.Abs(matrixH[lm1O + l].Real) < eps*tst1) + { + break; + } + + l--; + } + + var nm1 = n - 1; + var nm1O = nm1*order; + var nO = n*order; + var nOn = nO + n; + // Check for convergence + // One root found + if (l == n) + { + matrixH[nOn] += exshift; + vectorV[n] = matrixH[nOn]; + n--; + iter = 0; + } + else + { + // Form shift + Complex s; + if (iter != 10 && iter != 20) + { + s = matrixH[nOn]; + x = matrixH[nO + nm1]*matrixH[nm1O + n].Real; + + if (x.Real != 0.0 || x.Imaginary != 0.0) + { + y = (matrixH[nm1O + nm1] - 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[nm1O + n].Real) + Math.Abs(matrixH[(n - 2)*order + nm1].Real); + } + + for (var i = 0; i <= n; i++) + { + matrixH[i*order + i] -= s; + } + + exshift += s; + iter++; + + // Reduce to triangle (rows) + for (var i = l + 1; i <= n; i++) + { + var im1 = i - 1; + var im1O = im1*order; + var im1Oim1 = im1O + im1; + s = matrixH[im1O + i].Real; + norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real); + x = matrixH[im1Oim1]/norm; + vectorV[i - 1] = x; + matrixH[im1Oim1] = norm; + matrixH[im1O + i] = new Complex(0.0, s.Real/norm); + + for (var j = i; j < order; j++) + { + var jO = j*order; + y = matrixH[jO + im1]; + z = matrixH[jO + i]; + matrixH[jO + im1] = (x.Conjugate()*y) + (matrixH[im1O + i].Imaginary*z); + matrixH[jO + i] = (x*z) - (matrixH[im1O + i].Imaginary*y); + } + } + + s = matrixH[nOn]; + if (s.Imaginary != 0.0) + { + s /= matrixH[nOn].Magnitude; + matrixH[nOn] = matrixH[nOn].Magnitude; + + for (var j = n + 1; j < order; j++) + { + matrixH[j*order + n] *= s.Conjugate(); + } + } + + // Inverse operation (columns). + for (var j = l + 1; j <= n; j++) + { + x = vectorV[j - 1]; + var jO = j*order; + var jm1 = j - 1; + var jm1O = jm1*order; + var jm1Oj = jm1O + j; + for (var i = 0; i <= j; i++) + { + var jm1Oi = jm1O + i; + z = matrixH[jO + i]; + if (i != j) + { + y = matrixH[jm1Oi]; + matrixH[jm1Oi] = (x*y) + (matrixH[jm1O + j].Imaginary*z); + } + else + { + y = matrixH[jm1Oi].Real; + matrixH[jm1Oi] = new Complex((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[jm1O + j].Imaginary*z.Real), matrixH[jm1Oi].Imaginary); + } + + matrixH[jO + i] = (x.Conjugate()*z) - (matrixH[jm1O + j].Imaginary*y); + } + + for (var i = 0; i < order; i++) + { + y = dataEv[((j - 1)*order) + i]; + z = dataEv[(j*order) + i]; + dataEv[jm1O + i] = (x*y) + (matrixH[jm1Oj].Imaginary*z); + dataEv[jO + i] = (x.Conjugate()*z) - (matrixH[jm1Oj].Imaginary*y); + } + } + + if (s.Imaginary != 0.0) + { + for (var i = 0; i <= n; i++) + { + matrixH[nO + i] *= s; + } + + for (var i = 0; i < order; i++) + { + dataEv[nO + 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[j*order + i].Real) + Math.Abs(matrixH[j*order + i].Imaginary)); + } + } + + if (order == 1) + { + return; + } + + if (norm == 0.0) + { + return; + } + + for (n = order - 1; n > 0; n--) + { + var nO = n*order; + var nOn = nO + n; + x = vectorV[n]; + matrixH[nOn] = 1.0; + + for (var i = n - 1; i >= 0; i--) + { + z = 0.0; + for (var j = i + 1; j <= n; j++) + { + z += matrixH[j*order + i]*matrixH[nO + j]; + } + + y = x - vectorV[i]; + if (y.Real == 0.0 && y.Imaginary == 0.0) + { + y = eps*norm; + } + + matrixH[nO + i] = z/y; + + // Overflow control + var tr = Math.Abs(matrixH[nO + i].Real) + Math.Abs(matrixH[nO + i].Imaginary); + if ((eps*tr)*tr > 1) + { + for (var j = i; j <= n; j++) + { + matrixH[nO + j] = matrixH[nO + j]/tr; + } + } + } + } + + // Back transformation to get eigenvectors of original matrix + for (var j = order - 1; j > 0; j--) + { + var jO = j*order; + for (var i = 0; i < order; i++) + { + z = Complex.Zero; + for (var k = 0; k <= j; k++) + { + z += dataEv[(k*order) + i]*matrixH[jO + k]; + } + + dataEv[jO + i] = z; + } + } + } + /// /// Assumes that and have already been transposed. /// diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs index d9067ede..d2c8d351 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,12 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Complex32; -using MathNet.Numerics.LinearAlgebra.Complex32.Factorization; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.Managed { @@ -2487,9 +2485,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed var d = new float[order]; var e = new float[order]; - DenseEvd.SymmetricTridiagonalize(matrixCopy, d, e, tau, order); - DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); - DenseEvd.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); + SymmetricTridiagonalize(matrixCopy, d, e, tau, order); + SymmetricDiagonalize(matrixEv, d, e, order); + SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); for (var i = 0; i < order; i++) { @@ -2499,10 +2497,11 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed } else { - var v = new DenseVector(order); + var v = new Complex32[order]; + + NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); + NonsymmetricReduceHessenberToRealSchur(v, matrixEv, matrixCopy, order); - DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); - DenseEvd.NonsymmetricReduceHessenberToRealSchur(v.Values, matrixEv, matrixCopy, order); for (var i = 0; i < order; i++) { vectorEv[i] = new Complex(v[i].Real, v[i].Imaginary); @@ -2511,6 +2510,724 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed } } + /// + /// Reduces a complex Hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations. + /// + /// Source matrix to reduce + /// Output: Arrays for internal storage of real parts of eigenvalues + /// Output: Arrays for internal storage of imaginary parts of eigenvalues + /// Output: Arrays that contains further information about the transformations. + /// Order of initial matrix + /// 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. + internal 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*order + 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[k*order + i].Real) + Math.Abs(matrixA[k*order + i].Imaginary); + } + + if (scale == 0.0f) + { + tau[i - 1] = Complex32.One; + e[i] = 0.0f; + } + else + { + for (var k = 0; k < i; k++) + { + matrixA[k*order + i] /= scale; + h += matrixA[k*order + i].MagnitudeSquared; + } + + Complex32 g = (float) Math.Sqrt(h); + e[i] = scale*g.Real; + + Complex32 temp; + var im1Oi = (i - 1)*order + i; + var f = matrixA[im1Oi]; + if (f.Magnitude != 0.0f) + { + temp = -(matrixA[im1Oi].Conjugate()*tau[i].Conjugate())/f.Magnitude; + h += f.Magnitude*g.Real; + g = 1.0f + (g/f.Magnitude); + matrixA[im1Oi] *= g; + } + else + { + temp = -tau[i].Conjugate(); + matrixA[im1Oi] = g; + } + + if ((f.Magnitude == 0.0f) || (i != 1)) + { + f = Complex32.Zero; + for (var j = 0; j < i; j++) + { + var tmp = Complex32.Zero; + var jO = j*order; + // Form element of A*U. + for (var k = 0; k <= j; k++) + { + tmp += matrixA[k*order + j]*matrixA[k*order + i].Conjugate(); + } + + for (var k = j + 1; k <= i - 1; k++) + { + tmp += matrixA[jO + k].Conjugate()*matrixA[k*order + i].Conjugate(); + } + + // Form element of P + tau[j] = tmp/h; + f += (tmp/h)*matrixA[jO + i]; + } + + hh = f.Real/(h + h); + + // Form the reduced A. + for (var j = 0; j < i; j++) + { + f = matrixA[j*order + i].Conjugate(); + g = tau[j] - (hh*f); + tau[j] = g.Conjugate(); + + for (var k = 0; k <= j; k++) + { + matrixA[k*order + j] -= (f*tau[k]) + (g*matrixA[k*order + i]); + } + } + } + + for (var k = 0; k < i; k++) + { + matrixA[k*order + i] *= scale; + } + + tau[i - 1] = temp.Conjugate(); + } + + hh = d[i]; + d[i] = matrixA[i*order + i].Real; + matrixA[i*order + i] = new Complex32(hh, scale*(float) Math.Sqrt(h)); + } + + hh = d[0]; + d[0] = matrixA[0].Real; + matrixA[0] = hh; + e[0] = 0.0f; + } + + /// + /// Symmetric tridiagonal QL algorithm. + /// + /// Data array of matrix V (eigenvectors) + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + /// + internal 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.DoublePrecision; + 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 NonConvergenceException(); + } + } 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; + } + } + } + } + + /// + /// Determines eigenvectors by undoing the symmetric tridiagonalize transformation + /// + /// Data array of matrix V (eigenvectors) + /// Previously tridiagonalized matrix by SymmetricTridiagonalize. + /// Contains further information about the transformations + /// Input matrix order + /// 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. + internal 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*order + 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[k*order + i]; + } + + s = (s/h)/h; + + for (var k = 0; k < i; k++) + { + dataEv[(j*order) + k] -= s*matrixA[k*order + i].Conjugate(); + } + } + } + } + } + + /// + /// Nonsymmetric reduction to Hessenberg form. + /// + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Order of initial matrix + /// 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. + internal 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; + var mm1O = (m - 1)*order; + for (var i = m; i < order; i++) + { + scale += Math.Abs(matrixH[mm1O + i].Real) + Math.Abs(matrixH[mm1O + i].Imaginary); + } + + if (scale != 0.0f) + { + // Compute Householder transformation. + var h = 0.0f; + for (var i = order - 1; i >= m; i--) + { + ort[i] = matrixH[mm1O + i]/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[mm1O + m] = 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; + var jO = j*order; + for (var i = order - 1; i >= m; i--) + { + f += ort[i].Conjugate()*matrixH[jO + i]; + } + + f = f/h; + for (var i = m; i < order; i++) + { + matrixH[jO + i] -= 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[j*order + i]; + } + + f = f/h; + for (var j = m; j < order; j++) + { + matrixH[j*order + i] -= f*ort[j].Conjugate(); + } + } + + ort[m] = scale*ort[m]; + matrixH[mm1O + m] *= -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--) + { + var mm1O = (m - 1)*order; + var mm1Om = mm1O + m; + if (matrixH[mm1Om] != Complex32.Zero && ort[m] != Complex32.Zero) + { + var norm = (matrixH[mm1Om].Real*ort[m].Real) + (matrixH[mm1Om].Imaginary*ort[m].Imaginary); + + for (var i = m + 1; i < order; i++) + { + ort[i] = matrixH[mm1O + i]; + } + + 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++) + { + var im1 = i - 1; + var im1O = im1*order; + var im1Oi = im1O + i; + var iO = i*order; + if (matrixH[im1Oi].Imaginary != 0.0f) + { + var y = matrixH[im1Oi]/matrixH[im1Oi].Magnitude; + matrixH[im1Oi] = matrixH[im1Oi].Magnitude; + for (var j = i; j < order; j++) + { + matrixH[j*order + i] *= y.Conjugate(); + } + + for (var j = 0; j <= Math.Min(i + 1, order - 1); j++) + { + matrixH[iO + j] *= y; + } + + for (var j = 0; j < order; j++) + { + dataEv[(i*order) + j] *= y; + } + } + } + } + + /// + /// Nonsymmetric reduction from Hessenberg to real Schur form. + /// + /// Data array of the eigenvectors + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Order of initial matrix + /// 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. + internal static void NonsymmetricReduceHessenberToRealSchur(Complex32[] vectorV, Complex32[] dataEv, Complex32[] matrixH, int order) + { + // Initialize + var n = order - 1; + var eps = (float) Precision.SinglePrecision; + + float norm; + Complex32 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 lm1 = l - 1; + var lm1O = lm1*order; + var lO = l*order; + var tst1 = Math.Abs(matrixH[lm1O + lm1].Real) + Math.Abs(matrixH[lm1O + lm1].Imaginary) + Math.Abs(matrixH[lO + l].Real) + Math.Abs(matrixH[lO + l].Imaginary); + if (Math.Abs(matrixH[lm1O + l].Real) < eps*tst1) + { + break; + } + + l--; + } + + var nm1 = n - 1; + var nm1O = nm1*order; + var nO = n*order; + var nOn = nO + n; + // Check for convergence + // One root found + if (l == n) + { + matrixH[nOn] += exshift; + vectorV[n] = matrixH[nOn]; + n--; + iter = 0; + } + else + { + // Form shift + Complex32 s; + if (iter != 10 && iter != 20) + { + s = matrixH[nOn]; + x = matrixH[nO + nm1]*matrixH[nm1O + n].Real; + + if (x.Real != 0.0f || x.Imaginary != 0.0f) + { + y = (matrixH[nm1O + nm1] - s)/2.0f; + z = ((y*y) + x).SquareRoot(); + if ((y.Real*z.Real) + (y.Imaginary*z.Imaginary) < 0.0) + { + z *= -1.0f; + } + + x /= y + z; + s = s - x; + } + } + else + { + // Form exceptional shift + s = Math.Abs(matrixH[nm1O + n].Real) + Math.Abs(matrixH[(n - 2)*order + nm1].Real); + } + + for (var i = 0; i <= n; i++) + { + matrixH[i*order + i] -= s; + } + + exshift += s; + iter++; + + // Reduce to triangle (rows) + for (var i = l + 1; i <= n; i++) + { + var im1 = i - 1; + var im1O = im1*order; + var im1Oim1 = im1O + im1; + s = matrixH[im1O + i].Real; + norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real); + x = matrixH[im1Oim1]/norm; + vectorV[i - 1] = x; + matrixH[im1Oim1] = norm; + matrixH[im1O + i] = new Complex32(0.0f, s.Real/norm); + + for (var j = i; j < order; j++) + { + var jO = j*order; + y = matrixH[jO + im1]; + z = matrixH[jO + i]; + matrixH[jO + im1] = (x.Conjugate()*y) + (matrixH[im1O + i].Imaginary*z); + matrixH[jO + i] = (x*z) - (matrixH[im1O + i].Imaginary*y); + } + } + + s = matrixH[nOn]; + if (s.Imaginary != 0.0f) + { + s /= matrixH[nOn].Magnitude; + matrixH[nOn] = matrixH[nOn].Magnitude; + + for (var j = n + 1; j < order; j++) + { + matrixH[j*order + n] *= s.Conjugate(); + } + } + + // Inverse operation (columns). + for (var j = l + 1; j <= n; j++) + { + x = vectorV[j - 1]; + var jO = j*order; + var jm1 = j - 1; + var jm1O = jm1*order; + var jm1Oj = jm1O + j; + for (var i = 0; i <= j; i++) + { + var jm1Oi = jm1O + i; + z = matrixH[jO + i]; + if (i != j) + { + y = matrixH[jm1Oi]; + matrixH[jm1Oi] = (x*y) + (matrixH[jm1O + j].Imaginary*z); + } + else + { + y = matrixH[jm1Oi].Real; + matrixH[jm1Oi] = new Complex32((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[jm1O + j].Imaginary*z.Real), matrixH[jm1Oi].Imaginary); + } + + matrixH[jO + i] = (x.Conjugate()*z) - (matrixH[jm1O + j].Imaginary*y); + } + + for (var i = 0; i < order; i++) + { + y = dataEv[((j - 1)*order) + i]; + z = dataEv[(j*order) + i]; + dataEv[jm1O + i] = (x*y) + (matrixH[jm1Oj].Imaginary*z); + dataEv[jO + i] = (x.Conjugate()*z) - (matrixH[jm1Oj].Imaginary*y); + } + } + + if (s.Imaginary != 0.0f) + { + for (var i = 0; i <= n; i++) + { + matrixH[nO + i] *= s; + } + + for (var i = 0; i < order; i++) + { + dataEv[nO + 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[j*order + i].Real) + Math.Abs(matrixH[j*order + i].Imaginary)); + } + } + + if (order == 1) + { + return; + } + + if (norm == 0.0) + { + return; + } + + for (n = order - 1; n > 0; n--) + { + var nO = n*order; + var nOn = nO + n; + x = vectorV[n]; + matrixH[nOn] = 1.0f; + + for (var i = n - 1; i >= 0; i--) + { + z = 0.0f; + for (var j = i + 1; j <= n; j++) + { + z += matrixH[j*order + i]*matrixH[nO + j]; + } + + y = x - vectorV[i]; + if (y.Real == 0.0f && y.Imaginary == 0.0f) + { + y = eps*norm; + } + + matrixH[nO + i] = z/y; + + // Overflow control + var tr = Math.Abs(matrixH[nO + i].Real) + Math.Abs(matrixH[nO + i].Imaginary); + if ((eps*tr)*tr > 1) + { + for (var j = i; j <= n; j++) + { + matrixH[nO + j] = matrixH[nO + j]/tr; + } + } + } + } + + // Back transformation to get eigenvectors of original matrix + for (var j = order - 1; j > 0; j--) + { + var jO = j*order; + for (var i = 0; i < order; i++) + { + z = Complex32.Zero; + for (var k = 0; k <= j; k++) + { + z += dataEv[(k*order) + i]*matrixH[jO + k]; + } + + dataEv[jO + i] = z; + } + } + } + /// /// Assumes that and have already been transposed. /// diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs index 8ab8ee09..fa628798 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,10 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.Managed { @@ -2548,15 +2548,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed d[i] = matrixEv[i*order + om1]; } - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.SymmetricTridiagonalize(matrixEv, d, e, order); - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); + SymmetricTridiagonalize(matrixEv, d, e, order); + SymmetricDiagonalize(matrixEv, d, e, order); } else { var matrixH = new double[matrix.Length]; Buffer.BlockCopy(matrix, 0, matrixH, 0, matrix.Length*Constants.SizeOfDouble); - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); + NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); + NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); } for (var i = 0; i < order; i++) @@ -2577,5 +2577,1015 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed } } } + + /// + /// Symmetric Householder reduction to tridiagonal form. + /// + /// Data array of matrix V (eigenvectors) + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + internal static void SymmetricTridiagonalize(double[] a, double[] d, double[] e, int order) + { + // 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(d[k]); + } + + if (scale == 0.0) + { + e[i] = d[i - 1]; + for (var j = 0; j < i; j++) + { + d[j] = a[(j*order) + i - 1]; + a[(j*order) + i] = 0.0; + a[(i*order) + j] = 0.0; + } + } + else + { + // Generate Householder vector. + for (var k = 0; k < i; k++) + { + d[k] /= scale; + h += d[k]*d[k]; + } + + var f = d[i - 1]; + var g = Math.Sqrt(h); + if (f > 0) + { + g = -g; + } + + e[i] = scale*g; + h = h - (f*g); + d[i - 1] = f - g; + + for (var j = 0; j < i; j++) + { + e[j] = 0.0; + } + + // Apply similarity transformation to remaining columns. + for (var j = 0; j < i; j++) + { + f = d[j]; + a[(i*order) + j] = f; + g = e[j] + (a[(j*order) + j]*f); + + for (var k = j + 1; k <= i - 1; k++) + { + g += a[(j*order) + k]*d[k]; + e[k] += a[(j*order) + k]*f; + } + + e[j] = g; + } + + f = 0.0; + + for (var j = 0; j < i; j++) + { + e[j] /= h; + f += e[j]*d[j]; + } + + var hh = f/(h + h); + + for (var j = 0; j < i; j++) + { + e[j] -= hh*d[j]; + } + + for (var j = 0; j < i; j++) + { + f = d[j]; + g = e[j]; + + for (var k = j; k <= i - 1; k++) + { + a[(j*order) + k] -= (f*e[k]) + (g*d[k]); + } + + d[j] = a[(j*order) + i - 1]; + a[(j*order) + i] = 0.0; + } + } + + d[i] = h; + } + + // Accumulate transformations. + for (var i = 0; i < order - 1; i++) + { + a[(i*order) + order - 1] = a[(i*order) + i]; + a[(i*order) + i] = 1.0; + var h = d[i + 1]; + if (h != 0.0) + { + for (var k = 0; k <= i; k++) + { + d[k] = a[((i + 1)*order) + k]/h; + } + + for (var j = 0; j <= i; j++) + { + var g = 0.0; + for (var k = 0; k <= i; k++) + { + g += a[((i + 1)*order) + k]*a[(j*order) + k]; + } + + for (var k = 0; k <= i; k++) + { + a[(j*order) + k] -= g*d[k]; + } + } + } + + for (var k = 0; k <= i; k++) + { + a[((i + 1)*order) + k] = 0.0; + } + } + + for (var j = 0; j < order; j++) + { + d[j] = a[(j*order) + order - 1]; + a[(j*order) + order - 1] = 0.0; + } + + a[(order*order) - 1] = 1.0; + e[0] = 0.0; + } + + /// + /// Symmetric tridiagonal QL algorithm. + /// + /// Data array of matrix V (eigenvectors) + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + /// + internal static void SymmetricDiagonalize(double[] a, 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.DoublePrecision; + 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 = a[((i + 1)*order) + k]; + a[((i + 1)*order) + k] = (s*a[(i*order) + k]) + (c*h); + a[(i*order) + k] = (c*a[(i*order) + k]) - (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 NonConvergenceException(); + } + } 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 = a[(i*order) + j]; + a[(i*order) + j] = a[(k*order) + j]; + a[(k*order) + j] = p; + } + } + } + } + + /// + /// Nonsymmetric reduction to Hessenberg form. + /// + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Order of initial matrix + /// 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. + internal static void NonsymmetricReduceToHessenberg(double[] a, double[] matrixH, int order) + { + var ort = new double[order]; + var high = order - 1; + for (var m = 1; m <= high - 1; m++) + { + var mm1 = m - 1; + var mm1O = mm1*order; + // Scale column. + var scale = 0.0; + for (var i = m; i <= high; i++) + { + scale += Math.Abs(matrixH[mm1O + i]); + } + + if (scale != 0.0) + { + // Compute Householder transformation. + var h = 0.0; + for (var i = high; i >= m; i--) + { + ort[i] = matrixH[mm1O + i]/scale; + h += ort[i]*ort[i]; + } + + var g = 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 (var j = m; j < order; j++) + { + var jO = j*order; + var f = 0.0; + for (var i = order - 1; i >= m; i--) + { + f += ort[i]*matrixH[jO + i]; + } + + f = f/h; + + for (var i = m; i <= high; i++) + { + matrixH[jO + i] -= f*ort[i]; + } + } + + for (var i = 0; i <= high; i++) + { + var f = 0.0; + for (var j = high; j >= m; j--) + { + f += ort[j]*matrixH[j*order + i]; + } + f = f/h; + + for (var j = m; j <= high; j++) + { + matrixH[j*order + i] -= f*ort[j]; + } + } + + ort[m] = scale*ort[m]; + matrixH[mm1O + m] = scale*g; + } + } + + // Accumulate transformations (Algol's ortran). + for (var i = 0; i < order; i++) + { + for (var j = 0; j < order; j++) + { + a[(j*order) + i] = i == j ? 1.0 : 0.0; + } + } + + for (var m = high - 1; m >= 1; m--) + { + var mm1 = m - 1; + var mm1O = mm1*order; + var mm1Om = mm1O + m; + if (matrixH[mm1Om] != 0.0) + { + for (var i = m + 1; i <= high; i++) + { + ort[i] = matrixH[mm1O + i]; + } + + for (var j = m; j <= high; j++) + { + var g = 0.0; + var jO = j*order; + for (var i = m; i <= high; i++) + { + g += ort[i]*a[jO + i]; + } + + // Double division avoids possible underflow + g = (g/ort[m])/matrixH[mm1Om]; + + for (var i = m; i <= high; i++) + { + a[jO + i] += g*ort[i]; + } + } + } + } + } + + /// + /// Nonsymmetric reduction from Hessenberg to real Schur form. + /// + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + /// + internal static void NonsymmetricReduceHessenberToRealSchur(double[] a, double[] matrixH, double[] d, double[] e, int order) + { + // Initialize + var n = order - 1; + var eps = Math.Pow(2.0, -52.0); + var exshift = 0.0; + double p = 0, q = 0, r = 0, s = 0, z = 0; + double w, x, y; + + // Store roots isolated by balanc and compute matrix norm + var norm = 0.0; + for (var i = 0; i < order; i++) + { + for (var j = Math.Max(i - 1, 0); j < order; j++) + { + norm = norm + Math.Abs(matrixH[j*order + i]); + } + } + + // 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 lm1 = l - 1; + var lm1O = lm1*order; + s = Math.Abs(matrixH[lm1O + lm1]) + Math.Abs(matrixH[l*order + l]); + + if (s == 0.0) + { + s = norm; + } + + if (Math.Abs(matrixH[lm1O + l]) < eps*s) + { + break; + } + + l--; + } + + // Check for convergence + // One root found + if (l == n) + { + var index = n*order + n; + matrixH[index] += exshift; + d[n] = matrixH[index]; + e[n] = 0.0; + n--; + iter = 0; + + // Two roots found + } + else if (l == n - 1) + { + var nO = n*order; + var nm1 = n - 1; + var nm1O = nm1*order; + var nOn = nO + n; + + w = matrixH[nm1O + n]*matrixH[nO + nm1]; + p = (matrixH[nm1O + nm1] - matrixH[nOn])/2.0; + q = (p*p) + w; + z = Math.Sqrt(Math.Abs(q)); + + matrixH[nOn] += exshift; + matrixH[nm1O + nm1] += exshift; + x = matrixH[nOn]; + + // Real pair + if (q >= 0) + { + if (p >= 0) + { + z = p + z; + } + else + { + z = p - z; + } + + d[nm1] = x + z; + + d[n] = d[nm1]; + if (z != 0.0) + { + d[n] = x - (w/z); + } + + e[n - 1] = 0.0; + e[n] = 0.0; + x = matrixH[nm1O + n]; + s = Math.Abs(x) + Math.Abs(z); + p = x/s; + q = z/s; + r = Math.Sqrt((p*p) + (q*q)); + p = p/r; + q = q/r; + + // Row modification + for (var j = n - 1; j < order; j++) + { + var jO = j*order; + var jOn = jO + n; + z = matrixH[jO + nm1]; + matrixH[jO + nm1] = (q*z) + (p*matrixH[jOn]); + matrixH[jOn] = (q*matrixH[jOn]) - (p*z); + } + + // Column modification + for (var i = 0; i <= n; i++) + { + var nOi = nO + i; + z = matrixH[nm1O + i]; + matrixH[nm1O + i] = (q*z) + (p*matrixH[nOi]); + matrixH[nOi] = (q*matrixH[nOi]) - (p*z); + } + + // Accumulate transformations + for (var i = 0; i < order; i++) + { + var nOi = nO + i; + z = a[nm1O + i]; + a[nm1O + i] = (q*z) + (p*a[nOi]); + a[nOi] = (q*a[nOi]) - (p*z); + } + + // Complex pair + } + else + { + d[n - 1] = x + p; + d[n] = x + p; + e[n - 1] = z; + e[n] = -z; + } + + n = n - 2; + iter = 0; + + // No convergence yet + } + else + { + var nO = n*order; + var nm1 = n - 1; + var nm1O = nm1*order; + var nOn = nO + n; + + // Form shift + x = matrixH[nOn]; + y = 0.0; + w = 0.0; + if (l < n) + { + y = matrixH[nm1O + nm1]; + w = matrixH[nm1O + n]*matrixH[nO + nm1]; + } + + // Wilkinson's original ad hoc shift + if (iter == 10) + { + exshift += x; + for (var i = 0; i <= n; i++) + { + matrixH[i*order + i] -= x; + } + + s = Math.Abs(matrixH[nm1O + n]) + Math.Abs(matrixH[(n - 2)*order + nm1]); + x = y = 0.75*s; + w = (-0.4375)*s*s; + } + + // MATLAB's new ad hoc shift + if (iter == 30) + { + s = (y - x)/2.0; + s = (s*s) + w; + if (s > 0) + { + s = Math.Sqrt(s); + if (y < x) + { + s = -s; + } + + s = x - (w/(((y - x)/2.0) + s)); + for (var i = 0; i <= n; i++) + { + matrixH[i*order + i] -= s; + } + + exshift += s; + x = y = w = 0.964; + } + } + + iter = iter + 1; + if (iter >= 30*order) + { + throw new NonConvergenceException(); + } + + // Look for two consecutive small sub-diagonal elements + var m = n - 2; + while (m >= l) + { + var mp1 = m + 1; + var mm1 = m - 1; + var mO = m*order; + var mp1O = mp1*order; + var mm1O = mm1*order; + + z = matrixH[mO + m]; + r = x - z; + s = y - z; + p = (((r*s) - w)/matrixH[mO + mp1]) + matrixH[mp1O + m]; + q = matrixH[mp1O + mp1] - z - r - s; + r = matrixH[mp1O + (m + 2)]; + s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); + p = p/s; + q = q/s; + r = r/s; + + if (m == l) + { + break; + } + + if (Math.Abs(matrixH[mm1O + m])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[mm1O + mm1]) + Math.Abs(z) + Math.Abs(matrixH[mp1O + mp1])))) + { + break; + } + + m--; + } + + var mp2 = m + 2; + for (var i = mp2; i <= n; i++) + { + matrixH[(i - 2)*order + i] = 0.0; + if (i > mp2) + { + matrixH[(i - 3)*order + i] = 0.0; + } + } + + // Double QR step involving rows l:n and columns m:n + for (var k = m; k <= n - 1; k++) + { + var notlast = k != n - 1; + var kO = k*order; + var km1 = k - 1; + var kp1 = k + 1; + var kp2 = k + 2; + var kp1O = kp1*order; + var kp2O = kp2*order; + var km1O = km1*order; + if (k != m) + { + p = matrixH[km1O + k]; + q = matrixH[km1O + kp1]; + r = notlast ? matrixH[km1O + kp2] : 0.0; + x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); + if (x == 0.0) + { + continue; + } + + p = p/x; + q = q/x; + r = r/x; + } + + s = Math.Sqrt((p*p) + (q*q) + (r*r)); + + if (p < 0) + { + s = -s; + } + + if (s != 0.0) + { + if (k != m) + { + matrixH[km1O + k] = (-s)*x; + } + else if (l != m) + { + matrixH[km1O + k] = -matrixH[km1O + k]; + } + + p = p + s; + x = p/s; + y = q/s; + z = r/s; + q = q/p; + r = r/p; + + // Row modification + for (var j = k; j < order; j++) + { + var jO = j*order; + var jOk = jO + k; + var jOkp1 = jO + kp1; + var jOkp2 = jO + kp2; + p = matrixH[jOk] + (q*matrixH[jOkp1]); + if (notlast) + { + p = p + (r*matrixH[jOkp2]); + matrixH[jOkp2] -= (p*z); + } + + matrixH[jOk] -= (p*x); + matrixH[jOkp1] -= (p*y); + } + + // Column modification + for (var i = 0; i <= Math.Min(n, k + 3); i++) + { + p = (x*matrixH[kO + i]) + (y*matrixH[kp1O + i]); + + if (notlast) + { + p = p + (z*matrixH[kp2O + i]); + matrixH[kp2O + i] -= (p*r); + } + + matrixH[kO + i] -= p; + matrixH[kp1O + i] -= (p*q); + } + + // Accumulate transformations + for (var i = 0; i < order; i++) + { + p = (x*a[kO + i]) + (y*a[kp1O + i]); + + if (notlast) + { + p = p + (z*a[kp2O + i]); + a[kp2O + i] -= p*r; + } + + a[kO + i] -= p; + a[kp1O + i] -= p*q; + } + } // (s != 0) + } // k loop + } // check convergence + } // while (n >= low) + + // Backsubstitute to find vectors of upper triangular form + if (norm == 0.0) + { + return; + } + + for (n = order - 1; n >= 0; n--) + { + var nO = n*order; + var nm1 = n - 1; + var nm1O = nm1*order; + + p = d[n]; + q = e[n]; + + + // Real vector + double t; + if (q == 0.0) + { + var l = n; + matrixH[nO + n] = 1.0; + for (var i = n - 1; i >= 0; i--) + { + var ip1 = i + 1; + var iO = i*order; + var ip1O = ip1*order; + + w = matrixH[iO + i] - p; + r = 0.0; + for (var j = l; j <= n; j++) + { + r = r + (matrixH[j*order + i]*matrixH[nO + j]); + } + + if (e[i] < 0.0) + { + z = w; + s = r; + } + else + { + l = i; + if (e[i] == 0.0) + { + if (w != 0.0) + { + matrixH[nO + i] = (-r)/w; + } + else + { + matrixH[nO + i] = (-r)/(eps*norm); + } + + // Solve real equations + } + else + { + x = matrixH[ip1O + i]; + y = matrixH[iO + ip1]; + q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]); + t = ((x*s) - (z*r))/q; + matrixH[nO + i] = t; + if (Math.Abs(x) > Math.Abs(z)) + { + matrixH[nO + ip1] = (-r - (w*t))/x; + } + else + { + matrixH[nO + ip1] = (-s - (y*t))/z; + } + } + + // Overflow control + t = Math.Abs(matrixH[nO + i]); + if ((eps*t)*t > 1) + { + for (var j = i; j <= n; j++) + { + matrixH[nO + j] /= t; + } + } + } + } + + // Complex vector + } + else if (q < 0) + { + var l = n - 1; + + // Last vector component imaginary so matrix is triangular + if (Math.Abs(matrixH[nm1O + n]) > Math.Abs(matrixH[nO + nm1])) + { + matrixH[nm1O + nm1] = q/matrixH[nm1O + n]; + matrixH[nO + nm1] = (-(matrixH[nO + n] - p))/matrixH[nm1O + n]; + } + else + { + var res = Cdiv(0.0, -matrixH[nO + nm1], matrixH[nm1O + nm1] - p, q); + matrixH[nm1O + nm1] = res.Real; + matrixH[nO + nm1] = res.Imaginary; + } + + matrixH[nm1O + n] = 0.0; + matrixH[nO + n] = 1.0; + for (var i = n - 2; i >= 0; i--) + { + var ip1 = i + 1; + var iO = i*order; + var ip1O = ip1*order; + var ra = 0.0; + var sa = 0.0; + for (var j = l; j <= n; j++) + { + var jO = j*order; + var jOi = jO + i; + ra = ra + (matrixH[jOi]*matrixH[nm1O + j]); + sa = sa + (matrixH[jOi]*matrixH[nO + j]); + } + + w = matrixH[iO + i] - p; + + if (e[i] < 0.0) + { + z = w; + r = ra; + s = sa; + } + else + { + l = i; + if (e[i] == 0.0) + { + var res = Cdiv(-ra, -sa, w, q); + matrixH[nm1O + i] = res.Real; + matrixH[nO + i] = res.Imaginary; + } + else + { + // Solve complex equations + x = matrixH[ip1O + i]; + y = matrixH[iO + ip1]; + + var vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q); + var vi = (d[i] - p)*2.0*q; + if ((vr == 0.0) && (vi == 0.0)) + { + vr = eps*norm*(Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); + } + + var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi); + matrixH[nm1O + i] = res.Real; + matrixH[nO + i] = res.Imaginary; + if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q))) + { + matrixH[nm1O + ip1] = (-ra - (w*matrixH[nm1O + i]) + (q*matrixH[nO + i]))/x; + matrixH[nO + ip1] = (-sa - (w*matrixH[nO + i]) - (q*matrixH[nm1O + i]))/x; + } + else + { + res = Cdiv(-r - (y*matrixH[nm1O + i]), -s - (y*matrixH[nO + i]), z, q); + matrixH[nm1O + ip1] = res.Real; + matrixH[nO + ip1] = res.Imaginary; + } + } + + // Overflow control + t = Math.Max(Math.Abs(matrixH[nm1O + i]), Math.Abs(matrixH[nO + i])); + if ((eps*t)*t > 1) + { + for (var j = i; j <= n; j++) + { + matrixH[nm1O + j] /= t; + matrixH[nO + j] /= t; + } + } + } + } + } + } + + // Back transformation to get eigenvectors of original matrix + for (var j = order - 1; j >= 0; j--) + { + var jO = j*order; + for (var i = 0; i < order; i++) + { + z = 0.0; + for (var k = 0; k <= j; k++) + { + z = z + (a[k*order + i]*matrixH[jO + k]); + } + + a[jO + i] = z; + } + } + } + + /// + /// Complex scalar division X/Y. + /// + /// Real part of X + /// Imaginary part of X + /// Real part of Y + /// Imaginary part of Y + /// Division result as a number. + static Complex Cdiv(double xreal, double ximag, double yreal, double yimag) + { + if (Math.Abs(yimag) < Math.Abs(yreal)) + { + return new Complex((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal)))); + } + + return new Complex((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag)))); + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs index 9bb295f0..2258a484 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,10 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.Managed { @@ -2554,15 +2554,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed d[i] = matrixEv[i*order + om1]; } - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.SymmetricTridiagonalize(matrixEv, d, e, order); - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); + SymmetricTridiagonalize(matrixEv, d, e, order); + SymmetricDiagonalize(matrixEv, d, e, order); } else { var matrixH = new float[matrix.Length]; Buffer.BlockCopy(matrix, 0, matrixH, 0, matrix.Length*Constants.SizeOfFloat); - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); + NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); + NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); } for (var i = 0; i < order; i++) @@ -2582,5 +2582,1015 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed } } } + + /// + /// Symmetric Householder reduction to tridiagonal form. + /// + /// Data array of matrix V (eigenvectors) + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + internal static void SymmetricTridiagonalize(float[] a, float[] d, float[] e, int order) + { + // 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(d[k]); + } + + if (scale == 0.0f) + { + e[i] = d[i - 1]; + for (var j = 0; j < i; j++) + { + d[j] = a[(j*order) + i - 1]; + a[(j*order) + i] = 0.0f; + a[(i*order) + j] = 0.0f; + } + } + else + { + // Generate Householder vector. + for (var k = 0; k < i; k++) + { + d[k] /= scale; + h += d[k]*d[k]; + } + + var f = d[i - 1]; + var g = (float) Math.Sqrt(h); + if (f > 0) + { + g = -g; + } + + e[i] = scale*g; + h = h - (f*g); + d[i - 1] = f - g; + + for (var j = 0; j < i; j++) + { + e[j] = 0.0f; + } + + // Apply similarity transformation to remaining columns. + for (var j = 0; j < i; j++) + { + f = d[j]; + a[(i*order) + j] = f; + g = e[j] + (a[(j*order) + j]*f); + + for (var k = j + 1; k <= i - 1; k++) + { + g += a[(j*order) + k]*d[k]; + e[k] += a[(j*order) + k]*f; + } + + e[j] = g; + } + + f = 0.0f; + + for (var j = 0; j < i; j++) + { + e[j] /= h; + f += e[j]*d[j]; + } + + var hh = f/(h + h); + + for (var j = 0; j < i; j++) + { + e[j] -= hh*d[j]; + } + + for (var j = 0; j < i; j++) + { + f = d[j]; + g = e[j]; + + for (var k = j; k <= i - 1; k++) + { + a[(j*order) + k] -= (f*e[k]) + (g*d[k]); + } + + d[j] = a[(j*order) + i - 1]; + a[(j*order) + i] = 0.0f; + } + } + + d[i] = h; + } + + // Accumulate transformations. + for (var i = 0; i < order - 1; i++) + { + a[(i*order) + order - 1] = a[(i*order) + i]; + a[(i*order) + i] = 1.0f; + var h = d[i + 1]; + if (h != 0.0f) + { + for (var k = 0; k <= i; k++) + { + d[k] = a[((i + 1)*order) + k]/h; + } + + for (var j = 0; j <= i; j++) + { + var g = 0.0f; + for (var k = 0; k <= i; k++) + { + g += a[((i + 1)*order) + k]*a[(j*order) + k]; + } + + for (var k = 0; k <= i; k++) + { + a[(j*order) + k] -= g*d[k]; + } + } + } + + for (var k = 0; k <= i; k++) + { + a[((i + 1)*order) + k] = 0.0f; + } + } + + for (var j = 0; j < order; j++) + { + d[j] = a[(j*order) + order - 1]; + a[(j*order) + order - 1] = 0.0f; + } + + a[(order*order) - 1] = 1.0f; + e[0] = 0.0f; + } + + /// + /// Symmetric tridiagonal QL algorithm. + /// + /// Data array of matrix V (eigenvectors) + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + /// + internal static void SymmetricDiagonalize(float[] a, 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.SinglePrecision; + 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 = a[((i + 1)*order) + k]; + a[((i + 1)*order) + k] = (s*a[(i*order) + k]) + (c*h); + a[(i*order) + k] = (c*a[(i*order) + k]) - (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 NonConvergenceException(); + } + } 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 = a[(i*order) + j]; + a[(i*order) + j] = a[(k*order) + j]; + a[(k*order) + j] = p; + } + } + } + } + + /// + /// Nonsymmetric reduction to Hessenberg form. + /// + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Order of initial matrix + /// 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. + internal static void NonsymmetricReduceToHessenberg(float[] a, float[] matrixH, int order) + { + var ort = new float[order]; + var high = order - 1; + for (var m = 1; m <= high - 1; m++) + { + var mm1 = m - 1; + var mm1O = mm1*order; + // Scale column. + var scale = 0.0f; + for (var i = m; i <= high; i++) + { + scale += Math.Abs(matrixH[mm1O + i]); + } + + if (scale != 0.0f) + { + // Compute Householder transformation. + var h = 0.0f; + for (var i = high; i >= m; i--) + { + ort[i] = matrixH[mm1O + i]/scale; + h += ort[i]*ort[i]; + } + + var g = (float) 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 (var j = m; j < order; j++) + { + var jO = j*order; + var f = 0.0f; + for (var i = order - 1; i >= m; i--) + { + f += ort[i]*matrixH[jO + i]; + } + + f = f/h; + + for (var i = m; i <= high; i++) + { + matrixH[jO + i] -= f*ort[i]; + } + } + + for (var i = 0; i <= high; i++) + { + var f = 0.0f; + for (var j = high; j >= m; j--) + { + f += ort[j]*matrixH[j*order + i]; + } + f = f/h; + + for (var j = m; j <= high; j++) + { + matrixH[j*order + i] -= f*ort[j]; + } + } + + ort[m] = scale*ort[m]; + matrixH[mm1O + m] = scale*g; + } + } + + // Accumulate transformations (Algol's ortran). + for (var i = 0; i < order; i++) + { + for (var j = 0; j < order; j++) + { + a[(j*order) + i] = i == j ? 1.0f : 0.0f; + } + } + + for (var m = high - 1; m >= 1; m--) + { + var mm1 = m - 1; + var mm1O = mm1*order; + var mm1Om = mm1O + m; + if (matrixH[mm1Om] != 0.0) + { + for (var i = m + 1; i <= high; i++) + { + ort[i] = matrixH[mm1O + i]; + } + + for (var j = m; j <= high; j++) + { + var g = 0.0f; + var jO = j*order; + for (var i = m; i <= high; i++) + { + g += ort[i]*a[jO + i]; + } + + // Double division avoids possible underflow + g = (g/ort[m])/matrixH[mm1Om]; + + for (var i = m; i <= high; i++) + { + a[jO + i] += g*ort[i]; + } + } + } + } + } + + /// + /// Nonsymmetric reduction from Hessenberg to real Schur form. + /// + /// Data array of matrix V (eigenvectors) + /// Array for internal storage of nonsymmetric Hessenberg form. + /// Arrays for internal storage of real parts of eigenvalues + /// Arrays for internal storage of imaginary parts of eigenvalues + /// Order of initial matrix + /// 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. + /// + internal static void NonsymmetricReduceHessenberToRealSchur(float[] a, float[] matrixH, float[] d, float[] e, int order) + { + // Initialize + var n = order - 1; + var eps = (float) Precision.SinglePrecision; + var exshift = 0.0f; + float p = 0, q = 0, r = 0, s = 0, z = 0; + float w, x, y; + + // Store roots isolated by balanc and compute matrix norm + var norm = 0.0f; + for (var i = 0; i < order; i++) + { + for (var j = Math.Max(i - 1, 0); j < order; j++) + { + norm = norm + Math.Abs(matrixH[j*order + i]); + } + } + + // 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 lm1 = l - 1; + var lm1O = lm1*order; + s = Math.Abs(matrixH[lm1O + lm1]) + Math.Abs(matrixH[l*order + l]); + + if (s == 0.0) + { + s = norm; + } + + if (Math.Abs(matrixH[lm1O + l]) < eps*s) + { + break; + } + + l--; + } + + // Check for convergence + // One root found + if (l == n) + { + var index = n*order + n; + matrixH[index] += exshift; + d[n] = matrixH[index]; + e[n] = 0.0f; + n--; + iter = 0; + + // Two roots found + } + else if (l == n - 1) + { + var nO = n*order; + var nm1 = n - 1; + var nm1O = nm1*order; + var nOn = nO + n; + + w = matrixH[nm1O + n]*matrixH[nO + nm1]; + p = (matrixH[nm1O + nm1] - matrixH[nOn])/2.0f; + q = (p*p) + w; + z = (float) Math.Sqrt(Math.Abs(q)); + + matrixH[nOn] += exshift; + matrixH[nm1O + nm1] += exshift; + x = matrixH[nOn]; + + // Real pair + if (q >= 0) + { + if (p >= 0) + { + z = p + z; + } + else + { + z = p - z; + } + + d[nm1] = x + z; + + d[n] = d[nm1]; + if (z != 0.0) + { + d[n] = x - (w/z); + } + + e[n - 1] = 0.0f; + e[n] = 0.0f; + x = matrixH[nm1O + n]; + s = Math.Abs(x) + Math.Abs(z); + p = x/s; + q = z/s; + r = (float) Math.Sqrt((p*p) + (q*q)); + p = p/r; + q = q/r; + + // Row modification + for (var j = n - 1; j < order; j++) + { + var jO = j*order; + var jOn = jO + n; + z = matrixH[jO + nm1]; + matrixH[jO + nm1] = (q*z) + (p*matrixH[jOn]); + matrixH[jOn] = (q*matrixH[jOn]) - (p*z); + } + + // Column modification + for (var i = 0; i <= n; i++) + { + var nOi = nO + i; + z = matrixH[nm1O + i]; + matrixH[nm1O + i] = (q*z) + (p*matrixH[nOi]); + matrixH[nOi] = (q*matrixH[nOi]) - (p*z); + } + + // Accumulate transformations + for (var i = 0; i < order; i++) + { + var nOi = nO + i; + z = a[nm1O + i]; + a[nm1O + i] = (q*z) + (p*a[nOi]); + a[nOi] = (q*a[nOi]) - (p*z); + } + + // Complex pair + } + else + { + d[n - 1] = x + p; + d[n] = x + p; + e[n - 1] = z; + e[n] = -z; + } + + n = n - 2; + iter = 0; + + // No convergence yet + } + else + { + var nO = n*order; + var nm1 = n - 1; + var nm1O = nm1*order; + var nOn = nO + n; + + // Form shift + x = matrixH[nOn]; + y = 0.0f; + w = 0.0f; + if (l < n) + { + y = matrixH[nm1O + nm1]; + w = matrixH[nm1O + n]*matrixH[nO + nm1]; + } + + // Wilkinson's original ad hoc shift + if (iter == 10) + { + exshift += x; + for (var i = 0; i <= n; i++) + { + matrixH[i*order + i] -= x; + } + + s = Math.Abs(matrixH[nm1O + n]) + Math.Abs(matrixH[(n - 2)*order + nm1]); + x = y = 0.75f*s; + w = (-0.4375f)*s*s; + } + + // MATLAB's new ad hoc shift + if (iter == 30) + { + s = (y - x)/2.0f; + s = (s*s) + w; + if (s > 0) + { + s = (float) Math.Sqrt(s); + if (y < x) + { + s = -s; + } + + s = x - (w/(((y - x)/2.0f) + s)); + for (var i = 0; i <= n; i++) + { + matrixH[i*order + i] -= s; + } + + exshift += s; + x = y = w = 0.964f; + } + } + + iter = iter + 1; + if (iter >= 30*order) + { + throw new NonConvergenceException(); + } + + // Look for two consecutive small sub-diagonal elements + var m = n - 2; + while (m >= l) + { + var mp1 = m + 1; + var mm1 = m - 1; + var mO = m*order; + var mp1O = mp1*order; + var mm1O = mm1*order; + + z = matrixH[mO + m]; + r = x - z; + s = y - z; + p = (((r*s) - w)/matrixH[mO + mp1]) + matrixH[mp1O + m]; + q = matrixH[mp1O + mp1] - z - r - s; + r = matrixH[mp1O + (m + 2)]; + s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); + p = p/s; + q = q/s; + r = r/s; + + if (m == l) + { + break; + } + + if (Math.Abs(matrixH[mm1O + m])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[mm1O + mm1]) + Math.Abs(z) + Math.Abs(matrixH[mp1O + mp1])))) + { + break; + } + + m--; + } + + var mp2 = m + 2; + for (var i = mp2; i <= n; i++) + { + matrixH[(i - 2)*order + i] = 0.0f; + if (i > mp2) + { + matrixH[(i - 3)*order + i] = 0.0f; + } + } + + // Double QR step involving rows l:n and columns m:n + for (var k = m; k <= n - 1; k++) + { + var notlast = k != n - 1; + var kO = k*order; + var km1 = k - 1; + var kp1 = k + 1; + var kp2 = k + 2; + var kp1O = kp1*order; + var kp2O = kp2*order; + var km1O = km1*order; + if (k != m) + { + p = matrixH[km1O + k]; + q = matrixH[km1O + kp1]; + r = notlast ? matrixH[km1O + kp2] : 0.0f; + x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); + if (x == 0.0f) + { + continue; + } + + p = p/x; + q = q/x; + r = r/x; + } + + s = (float) Math.Sqrt((p*p) + (q*q) + (r*r)); + + if (p < 0) + { + s = -s; + } + + if (s != 0.0f) + { + if (k != m) + { + matrixH[km1O + k] = (-s)*x; + } + else if (l != m) + { + matrixH[km1O + k] = -matrixH[km1O + k]; + } + + p = p + s; + x = p/s; + y = q/s; + z = r/s; + q = q/p; + r = r/p; + + // Row modification + for (var j = k; j < order; j++) + { + var jO = j*order; + var jOk = jO + k; + var jOkp1 = jO + kp1; + var jOkp2 = jO + kp2; + p = matrixH[jOk] + (q*matrixH[jOkp1]); + if (notlast) + { + p = p + (r*matrixH[jOkp2]); + matrixH[jOkp2] -= (p*z); + } + + matrixH[jOk] -= (p*x); + matrixH[jOkp1] -= (p*y); + } + + // Column modification + for (var i = 0; i <= Math.Min(n, k + 3); i++) + { + p = (x*matrixH[kO + i]) + (y*matrixH[kp1O + i]); + + if (notlast) + { + p = p + (z*matrixH[kp2O + i]); + matrixH[kp2O + i] -= (p*r); + } + + matrixH[kO + i] -= p; + matrixH[kp1O + i] -= (p*q); + } + + // Accumulate transformations + for (var i = 0; i < order; i++) + { + p = (x*a[kO + i]) + (y*a[kp1O + i]); + + if (notlast) + { + p = p + (z*a[kp2O + i]); + a[kp2O + i] -= p*r; + } + + a[kO + i] -= p; + a[kp1O + i] -= p*q; + } + } // (s != 0) + } // k loop + } // check convergence + } // while (n >= low) + + // Backsubstitute to find vectors of upper triangular form + if (norm == 0.0f) + { + return; + } + + for (n = order - 1; n >= 0; n--) + { + var nO = n*order; + var nm1 = n - 1; + var nm1O = nm1*order; + + p = d[n]; + q = e[n]; + + + // Real vector + float t; + if (q == 0.0f) + { + var l = n; + matrixH[nO + n] = 1.0f; + for (var i = n - 1; i >= 0; i--) + { + var ip1 = i + 1; + var iO = i*order; + var ip1O = ip1*order; + + w = matrixH[iO + i] - p; + r = 0.0f; + for (var j = l; j <= n; j++) + { + r = r + (matrixH[j*order + i]*matrixH[nO + j]); + } + + if (e[i] < 0.0) + { + z = w; + s = r; + } + else + { + l = i; + if (e[i] == 0.0f) + { + if (w != 0.0f) + { + matrixH[nO + i] = (-r)/w; + } + else + { + matrixH[nO + i] = (-r)/(eps*norm); + } + + // Solve real equations + } + else + { + x = matrixH[ip1O + i]; + y = matrixH[iO + ip1]; + q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]); + t = ((x*s) - (z*r))/q; + matrixH[nO + i] = t; + if (Math.Abs(x) > Math.Abs(z)) + { + matrixH[nO + ip1] = (-r - (w*t))/x; + } + else + { + matrixH[nO + ip1] = (-s - (y*t))/z; + } + } + + // Overflow control + t = Math.Abs(matrixH[nO + i]); + if ((eps*t)*t > 1) + { + for (var j = i; j <= n; j++) + { + matrixH[nO + j] /= t; + } + } + } + } + + // Complex vector + } + else if (q < 0) + { + var l = n - 1; + + // Last vector component imaginary so matrix is triangular + if (Math.Abs(matrixH[nm1O + n]) > Math.Abs(matrixH[nO + nm1])) + { + matrixH[nm1O + nm1] = q/matrixH[nm1O + n]; + matrixH[nO + nm1] = (-(matrixH[nO + n] - p))/matrixH[nm1O + n]; + } + else + { + var res = Cdiv(0.0f, -matrixH[nO + nm1], matrixH[nm1O + nm1] - p, q); + matrixH[nm1O + nm1] = res.Real; + matrixH[nO + nm1] = res.Imaginary; + } + + matrixH[nm1O + n] = 0.0f; + matrixH[nO + n] = 1.0f; + for (var i = n - 2; i >= 0; i--) + { + var ip1 = i + 1; + var iO = i*order; + var ip1O = ip1*order; + var ra = 0.0f; + var sa = 0.0f; + for (var j = l; j <= n; j++) + { + var jO = j*order; + var jOi = jO + i; + ra = ra + (matrixH[jOi]*matrixH[nm1O + j]); + sa = sa + (matrixH[jOi]*matrixH[nO + j]); + } + + w = matrixH[iO + i] - p; + + if (e[i] < 0.0) + { + z = w; + r = ra; + s = sa; + } + else + { + l = i; + if (e[i] == 0.0) + { + var res = Cdiv(-ra, -sa, w, q); + matrixH[nm1O + i] = res.Real; + matrixH[nO + i] = res.Imaginary; + } + else + { + // Solve complex equations + x = matrixH[ip1O + i]; + y = matrixH[iO + ip1]; + + var vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q); + var vi = (d[i] - p)*2.0f*q; + if ((vr == 0.0f) && (vi == 0.0f)) + { + vr = eps*norm*(Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); + } + + var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi); + matrixH[nm1O + i] = res.Real; + matrixH[nO + i] = res.Imaginary; + if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q))) + { + matrixH[nm1O + ip1] = (-ra - (w*matrixH[nm1O + i]) + (q*matrixH[nO + i]))/x; + matrixH[nO + ip1] = (-sa - (w*matrixH[nO + i]) - (q*matrixH[nm1O + i]))/x; + } + else + { + res = Cdiv(-r - (y*matrixH[nm1O + i]), -s - (y*matrixH[nO + i]), z, q); + matrixH[nm1O + ip1] = res.Real; + matrixH[nO + ip1] = res.Imaginary; + } + } + + // Overflow control + t = Math.Max(Math.Abs(matrixH[nm1O + i]), Math.Abs(matrixH[nO + i])); + if ((eps*t)*t > 1) + { + for (var j = i; j <= n; j++) + { + matrixH[nm1O + j] /= t; + matrixH[nO + j] /= t; + } + } + } + } + } + } + + // Back transformation to get eigenvectors of original matrix + for (var j = order - 1; j >= 0; j--) + { + var jO = j*order; + for (var i = 0; i < order; i++) + { + z = 0.0f; + for (var k = 0; k <= j; k++) + { + z = z + (a[k*order + i]*matrixH[jO + k]); + } + + a[jO + i] = z; + } + } + } + + /// + /// Complex scalar division X/Y. + /// + /// Real part of X + /// Imaginary part of X + /// Real part of Y + /// Imaginary part of Y + /// Division result as a number. + static Complex32 Cdiv(float xreal, float ximag, float yreal, float yimag) + { + if (Math.Abs(yimag) < Math.Abs(yreal)) + { + return new Complex32((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal)))); + } + + return new Complex32((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag)))); + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs index 93e5e7df..8c89a6b6 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -68,7 +68,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed /// /// Assumes that and have already been transposed. /// - protected static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, T[] matrix, T[] row) + static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, T[] matrix, T[] row) { if (transpose == Transpose.DontTranspose) { @@ -86,7 +86,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed /// /// Assumes that and have already been transposed. /// - protected static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, T[] matrix, T[] column) + static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, T[] matrix, T[] column) { if (transpose == Transpose.DontTranspose) { diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs index d01e94f1..d5380a6a 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,11 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Complex.Factorization; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference { @@ -2743,9 +2742,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference var d = new double[order]; var e = new double[order]; - DenseEvd.SymmetricTridiagonalize(matrixCopy, d, e, tau, order); - DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); - DenseEvd.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); + Managed.ManagedLinearAlgebraProvider.SymmetricTridiagonalize(matrixCopy, d, e, tau, order); + Managed.ManagedLinearAlgebraProvider.SymmetricDiagonalize(matrixEv, d, e, order); + Managed.ManagedLinearAlgebraProvider.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); for (var i = 0; i < order; i++) { @@ -2754,8 +2753,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference } else { - DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); - DenseEvd.NonsymmetricReduceHessenberToRealSchur(vectorEv, matrixEv, matrixCopy, order); + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceHessenberToRealSchur(vectorEv, matrixEv, matrixCopy, order); } for (var i = 0; i < order; i++) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs index 1a4575fa..ca19ab35 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,12 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Complex32; -using MathNet.Numerics.LinearAlgebra.Complex32.Factorization; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference { @@ -2743,9 +2741,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference var d = new float[order]; var e = new float[order]; - DenseEvd.SymmetricTridiagonalize(matrixCopy, d, e, tau, order); - DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); - DenseEvd.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); + Managed.ManagedLinearAlgebraProvider.SymmetricTridiagonalize(matrixCopy, d, e, tau, order); + Managed.ManagedLinearAlgebraProvider.SymmetricDiagonalize(matrixEv, d, e, order); + Managed.ManagedLinearAlgebraProvider.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order); for (var i = 0; i < order; i++) { @@ -2755,10 +2753,11 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference } else { - var v = new DenseVector(order); + var v = new Complex32[order]; + + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceHessenberToRealSchur(v, matrixEv, matrixCopy, order); - DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order); - DenseEvd.NonsymmetricReduceHessenberToRealSchur(v.Values, matrixEv, matrixCopy, order); for (var i = 0; i < order; i++) { vectorEv[i] = new Complex(v[i].Real, v[i].Imaginary); diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs index 3b12c727..a7969245 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,10 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference { @@ -2692,15 +2692,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference d[i] = matrixEv[i*order + om1]; } - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.SymmetricTridiagonalize(matrixEv, d, e, order); - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); + Managed.ManagedLinearAlgebraProvider.SymmetricTridiagonalize(matrixEv, d, e, order); + Managed.ManagedLinearAlgebraProvider.SymmetricDiagonalize(matrixEv, d, e, order); } else { var matrixH = new double[matrix.Length]; Buffer.BlockCopy(matrix, 0, matrixH, 0, matrix.Length*Constants.SizeOfDouble); - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); - Numerics.LinearAlgebra.Double.Factorization.DenseEvd.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); } for (var i = 0; i < order; i++) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs index f0dda826..c21e7520 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,10 +28,10 @@ // using System; -using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; +using QRMethod = MathNet.Numerics.LinearAlgebra.Factorization.QRMethod; namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference { @@ -2698,15 +2698,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference d[i] = matrixEv[i*order + om1]; } - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.SymmetricTridiagonalize(matrixEv, d, e, order); - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order); + Managed.ManagedLinearAlgebraProvider.SymmetricTridiagonalize(matrixEv, d, e, order); + Managed.ManagedLinearAlgebraProvider.SymmetricDiagonalize(matrixEv, d, e, order); } else { var matrixH = new float[matrix.Length]; Buffer.BlockCopy(matrix, 0, matrixH, 0, matrix.Length*Constants.SizeOfFloat); - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); - Numerics.LinearAlgebra.Single.Factorization.DenseEvd.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order); + Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order); } for (var i = 0; i < order; i++) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs index 2cb29b40..3fa70a14 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2020 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -27,8 +27,6 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; - namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference { /// @@ -64,41 +62,5 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference { return "Managed"; } - - /// - /// Assumes that and have already been transposed. - /// - protected static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, T[] matrix, T[] row) - { - if (transpose == Transpose.DontTranspose) - { - for (int i = 0; i < numCols; i++) - { - row[i] = matrix[(i * numRows) + rowindx]; - } - } - else - { - Array.Copy(matrix, rowindx * numCols, row, 0, numCols); - } - } - - /// - /// Assumes that and have already been transposed. - /// - protected static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, T[] matrix, T[] column) - { - if (transpose == Transpose.DontTranspose) - { - Array.Copy(matrix, colindx * numRows, column, 0, numRows); - } - else - { - for (int i = 0; i < numRows; i++) - { - column[i] = matrix[(i * numCols) + colindx]; - } - } - } } }