using System; namespace Science.Mathematics.LinearAlgebra { public class FactorizationAeqQLQt { private Matrix Q; private Matrix L; private double[,] a; public FactorizationAeqQLQt(SymmetricMatrix A) { a = new double[A.SizeOfRow, A.SizeOfColumn]; for (int i = 0; i < a.GetLength(0); i++) for (int j = 0; j < a.GetLength(1); j++) a[i, j] = A[i, j]; Q = new Matrix(A.SizeOfRow,A.SizeOfColumn); L = new Matrix(A.SizeOfRow,A.SizeOfColumn); } public void Compute() { Symmeig obj = new Symmeig(a); for (int i = 0; i < a.GetLength(0); i++) { L[i, i] = obj.d[i]; for (int j = 0; j < a.GetLength(1); j++) { Q[i, j] = obj.z[i, j]; } } } public Matrix Eigenvectors { get { return Q; } } public Matrix Eigenvalues { get { return L; } } } class Jacobi { public void eigsrt(double[] d) { eigsrt(d, null); } public void eigsrt(double[] d, double[,] v) { int k; int n = d.Length; for (int i = 0; i < n - 1; i++) { double p = d[k = i]; for (int j = i; j < n; j++) if (d[j] >= p) p = d[k = j]; if (k != i) { d[k] = d[i]; d[i] = p; if (v != null) for (int j = 0; j < n; j++) { p = v[j, i]; v[j, i] = v[j, k]; v[j, k] = p; } } } } int n; double[,] a; public double[,] v; public double[] d; public int nrot; double EPS; public Jacobi() { } public Jacobi(double[,] aa) { n = aa.GetLength(0); a = aa; v = new double[n, n]; d = new double[n]; nrot = 0; EPS = 1.0E-16; int i, j, ip, iq; double tresh, theta, tau, t, sm, s, h, g, c; double[] b = new double[n]; double[] z = new double[n]; for (ip = 0; ip < n; ip++) { for (iq = 0; iq < n; iq++) v[ip, iq] = 0.0; v[ip, ip] = 1.0; } for (ip = 0; ip < n; ip++) { b[ip] = d[ip] = a[ip, ip]; z[ip] = 0.0; } for (i = 1; i <= 50; i++) { sm = 0.0; for (ip = 0; ip < n - 1; ip++) { for (iq = ip + 1; iq < n; iq++) sm += Math.Abs(a[ip, iq]); } if (sm == 0.0) { eigsrt(d, v); return; } if (i < 4) tresh = 0.2 * sm / (n * n); else tresh = 0.0; for (ip = 0; ip < n - 1; ip++) { for (iq = ip + 1; iq < n; iq++) { g = 100.0 * Math.Abs(a[ip, iq]); if (i > 4 && g <= EPS * Math.Abs(d[ip]) && g <= EPS * Math.Abs(d[iq])) a[ip, iq] = 0.0; else if (Math.Abs(a[ip, iq]) > tresh) { h = d[iq] - d[ip]; if (g <= EPS * Math.Abs(h)) t = (a[ip, iq]) / h; else { theta = 0.5 * h / (a[ip, iq]); t = 1.0 / (Math.Abs(theta) + Math.Sqrt(1.0 + theta * theta)); if (theta < 0.0) t = -t; } c = 1.0 / Math.Sqrt(1 + t * t); s = t * c; tau = s / (1.0 + c); h = t * a[ip, iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip, iq] = 0.0; for (j = 0; j < ip; j++) rot(a, s, tau, j, ip, j, iq); for (j = ip + 1; j < iq; j++) rot(a, s, tau, ip, j, j, iq); for (j = iq + 1; j < n; j++) rot(a, s, tau, ip, j, iq, j); for (j = 0; j < n; j++) rot(v, s, tau, j, ip, j, iq); ++nrot; } } } for (ip = 0; ip < n; ip++) { b[ip] += z[ip]; d[ip] = b[ip]; z[ip] = 0.0; } } try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Too many iterations in routine jacobi"); } } private void rot(double[,] a, double s, double tau, int i, int j, int k, int l) { double g = a[i, j]; double h = a[k, l]; a[i, j] = g - s * (h + g * tau); a[k, l] = h + s * (g - h * tau); } } class Symmeig { int n; public double[,] z; public double[] d; double[] e; bool yesvecs; public Symmeig(double[,] a) { MSymmeig(a, true); } public Symmeig(double[,] a, bool yesvec) { MSymmeig(a, yesvec); } private void MSymmeig(double[,] a, bool yesvec) { n = a.GetLength(0); z = a; d = new double[n]; e = new double[n]; yesvecs = yesvec; tred2(); tqli(); sort(); } public Symmeig(double[] dd, double[] ee) { MMSymmeig(dd, ee, true); } public Symmeig(double[] dd, double[] ee, bool yesvec) { MMSymmeig(dd, ee, yesvec); } private void MMSymmeig(double[] dd, double[] ee, bool yesvec) { n = dd.Length; d = dd; e = ee; z = new double[n, n]; yesvecs = yesvec; for (int i = 0; i < n; i++) z[i, i] = 1.0; tqli(); sort(); } private void sort() { Jacobi obj = new Jacobi(); if (yesvecs) obj.eigsrt(d, z); else obj.eigsrt(d); } private void tred2() { int l, k, j, i; double scale, hh, h, g, f; for (i = n - 1; i > 0; i--) { l = i - 1; h = scale = 0.0; if (l > 0) { for (k = 0; k < i; k++) scale += Math.Abs(z[i, k]); if (scale == 0.0) e[i] = z[i, l]; else { for (k = 0; k < i; k++) { z[i, k] /= scale; h += z[i, k] * z[i, k]; } f = z[i, l]; g = (f >= 0.0 ? -Math.Sqrt(h) : Math.Sqrt(h)); e[i] = scale * g; h -= f * g; z[i, l] = f - g; f = 0.0; for (j = 0; j < i; j++) { if (yesvecs) z[j, i] = z[i, j] / h; g = 0.0; for (k = 0; k < j + 1; k++) g += z[j, k] * z[i, k]; for (k = j + 1; k < i; k++) g += z[k, j] * z[i, k]; e[j] = g / h; f += e[j] * z[i, j]; } hh = f / (h + h); for (j = 0; j < i; j++) { f = z[i, j]; e[j] = g = e[j] - hh * f; for (k = 0; k < j + 1; k++) z[j, k] -= (f * e[k] + g * z[i, k]); } } } else e[i] = z[i, l]; d[i] = h; } if (yesvecs) d[0] = 0.0; e[0] = 0.0; for (i = 0; i < n; i++) { if (yesvecs) { if (d[i] != 0.0) { for (j = 0; j < i; j++) { g = 0.0; for (k = 0; k < i; k++) g += z[i, k] * z[k, j]; for (k = 0; k < i; k++) z[k, j] -= g * z[k, i]; } } d[i] = z[i, i]; z[i, i] = 1.0; for (j = 0; j < i; j++) z[j, i] = z[i, j] = 0.0; } else { d[i] = z[i, i]; } } } private void tqli() { int m, l, iter, i, k; double s, r, p, g, f, dd, c, b; double EPS = 1.0E-16; for (i = 1; i < n; i++) e[i - 1] = e[i]; e[n - 1] = 0.0; for (l = 0; l < n; l++) { iter = 0; do { for (m = l; m < n - 1; m++) { dd = Math.Abs(d[m]) + Math.Abs(d[m + 1]); if (Math.Abs(e[m]) <= EPS * dd) break; } if (m != l) { if (iter++ == 30) try { throw new Science.Error(); } catch (Science.Error ee) { ee.Write("Too many iterations in tqli"); } g = (d[l + 1] - d[l]) / (2.0 * e[l]); r = pythag(g, 1.0); g = d[m] - d[l] + e[l] / (g + SIGN(r, g)); s = c = 1.0; p = 0.0; for (i = m - 1; i >= l; i--) { f = s * e[i]; b = c * e[i]; e[i + 1] = (r = pythag(f, g)); if (r == 0.0) { d[i + 1] -= p; e[m] = 0.0; break; } s = f / r; c = g / r; g = d[i + 1] - p; r = (d[i] - g) * s + 2.0 * c * b; d[i + 1] = g + (p = s * r); g = c * r - b; if (yesvecs) { for (k = 0; k < n; k++) { f = z[k, i + 1]; z[k, i + 1] = s * z[k, i] + c * f; z[k, i] = c * z[k, i] - s * f; } } } if (r == 0.0 && i >= l) continue; d[l] -= p; e[l] = g; e[m] = 0.0; } } while (m != l); } } private double pythag(double a, double b) { double absa = Math.Abs(a), absb = Math.Abs(b); return (absa > absb ? absa * Math.Sqrt(1.0 + SQR(absb / absa)) : (absb == 0.0 ? 0.0 : absb * Math.Sqrt(1.0 + SQR(absa / absb)))); } private static double SQR(double a) { return a * a; } private static double SIGN(double a, double b) { return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); } } }