Created
November 8, 2014 04:53
-
-
Save lh3/c280b2ac477e85c5c666 to your computer and use it in GitHub Desktop.
Eigenvalues/vectors for symmetric and asymmetric matrices
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #include <math.h> | |
| #include <stdio.h> | |
| #include <stdlib.h> | |
| #define RADIX 2.0 | |
| /************************* | |
| * balance a real matrix * | |
| *************************/ | |
| static void balanc(double **a, int n) | |
| { | |
| int i, j, last = 0; | |
| double s, r, g, f, c, sqrdx; | |
| sqrdx = RADIX * RADIX; | |
| while (last == 0) { | |
| last = 1; | |
| for (i = 0; i < n; i++) { | |
| r = c = 0.0; | |
| for (j = 0; j < n; j++) | |
| if (j != i) { | |
| c += fabs(a[j][i]); | |
| r += fabs(a[i][j]); | |
| } | |
| if (c != 0.0 && r != 0.0) { | |
| g = r / RADIX; | |
| f = 1.0; | |
| s = c + r; | |
| while (c < g) { | |
| f *= RADIX; | |
| c *= sqrdx; | |
| } | |
| g = r * RADIX; | |
| while (c > g) { | |
| f /= RADIX; | |
| c /= sqrdx; | |
| } | |
| if ((c + r) / f < 0.95 * s) { | |
| last = 0; | |
| g = 1.0 / f; | |
| for (j = 0; j < n; j++) | |
| a[i][j] *= g; | |
| for (j = 0; j < n; j++) | |
| a[j][i] *= f; | |
| } | |
| } | |
| } | |
| } | |
| } | |
| #define SWAP(a,b) do { double t = (a); (a) = (b); (b) = t; } while (0) | |
| /***************************************************** | |
| * convert a non-symmetric matrix to Hessenberg form * | |
| *****************************************************/ | |
| static void elmhes(double **a, int n) | |
| { | |
| int i, j, m; | |
| double y, x; | |
| for (m = 1; m < n - 1; m++) { | |
| x = 0.0; | |
| i = m; | |
| for (j = m; j < n; j++) { | |
| if (fabs(a[j][m - 1]) > fabs(x)) { | |
| x = a[j][m - 1]; | |
| i = j; | |
| } | |
| } | |
| if (i != m) { | |
| for (j = m - 1; j < n; j++) | |
| SWAP(a[i][j], a[m][j]); | |
| for (j = 0; j < n; j++) | |
| SWAP(a[j][i], a[j][m]); | |
| } | |
| if (x != 0.0) { | |
| for (i = m + 1; i < n; i++) { | |
| if ((y = a[i][m - 1]) != 0.0) { | |
| y /= x; | |
| a[i][m - 1] = y; | |
| for (j = m; j < n; j++) | |
| a[i][j] -= y * a[m][j]; | |
| for (j = 0; j < n; j++) | |
| a[j][m] += y * a[j][i]; | |
| } | |
| } | |
| } | |
| } | |
| } | |
| #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) | |
| /************************************** | |
| * QR algorithm for Hessenberg matrix * | |
| **************************************/ | |
| static void hqr(double **a, int n, double *wr, double *wi) | |
| { | |
| int nn, m, l, k, j, its, i, mmin; | |
| double z, y, x, w, v, u, t, s, r, q, p, anorm; | |
| p = q = r = 0.0; | |
| anorm = 0.0; | |
| for (i = 0; i < n; i++) | |
| for (j = i - 1 > 0 ? i - 1 : 0; j < n; j++) | |
| anorm += fabs(a[i][j]); | |
| nn = n - 1; | |
| t = 0.0; | |
| while (nn >= 0) { | |
| its = 0; | |
| do { | |
| for (l = nn; l > 0; l--) { | |
| s = fabs(a[l - 1][l - 1]) + fabs(a[l][l]); | |
| if (s == 0.0) | |
| s = anorm; | |
| if (fabs(a[l][l - 1]) + s == s) { | |
| a[l][l - 1] = 0.0; | |
| break; | |
| } | |
| } | |
| x = a[nn][nn]; | |
| if (l == nn) { | |
| wr[nn] = x + t; | |
| wi[nn--] = 0.0; | |
| } else { | |
| y = a[nn - 1][nn - 1]; | |
| w = a[nn][nn - 1] * a[nn - 1][nn]; | |
| if (l == nn - 1) { | |
| p = 0.5 * (y - x); | |
| q = p * p + w; | |
| z = sqrt(fabs(q)); | |
| x += t; | |
| if (q >= 0.0) { | |
| z = p + SIGN(z, p); | |
| wr[nn - 1] = wr[nn] = x + z; | |
| if (z != 0.0) | |
| wr[nn] = x - w / z; | |
| wi[nn - 1] = wi[nn] = 0.0; | |
| } else { | |
| wr[nn - 1] = wr[nn] = x + p; | |
| wi[nn - 1] = -(wi[nn] = z); | |
| } | |
| nn -= 2; | |
| } else { | |
| if (its == 30) { | |
| fprintf(stderr, "[hqr] too many iterations.\n"); | |
| break; | |
| } | |
| if (its == 10 || its == 20) { | |
| t += x; | |
| for (i = 0; i < nn + 1; i++) | |
| a[i][i] -= x; | |
| s = fabs(a[nn][nn - 1]) + fabs(a[nn - 1][nn - 2]); | |
| y = x = 0.75 * s; | |
| w = -0.4375 * s * s; | |
| } | |
| ++its; | |
| for (m = nn - 2; m >= l; m--) { | |
| z = a[m][m]; | |
| r = x - z; | |
| s = y - z; | |
| p = (r * s - w) / a[m + 1][m] + a[m][m + 1]; | |
| q = a[m + 1][m + 1] - z - r - s; | |
| r = a[m + 2][m + 1]; | |
| s = fabs(p) + fabs(q) + fabs(r); | |
| p /= s; | |
| q /= s; | |
| r /= s; | |
| if (m == l) | |
| break; | |
| u = fabs(a[m][m - 1]) * (fabs(q) + fabs(r)); | |
| v = fabs(p) * (fabs(a[m - 1][m - 1]) + fabs(z) + fabs(a[m + 1][m + 1])); | |
| if (u + v == v) | |
| break; | |
| } | |
| for (i = m; i < nn - 1; i++) { | |
| a[i + 2][i] = 0.0; | |
| if (i != m) | |
| a[i + 2][i - 1] = 0.0; | |
| } | |
| for (k = m; k < nn; k++) { | |
| if (k != m) { | |
| p = a[k][k - 1]; | |
| q = a[k + 1][k - 1]; | |
| r = 0.0; | |
| if (k + 1 != nn) | |
| r = a[k + 2][k - 1]; | |
| if ((x = fabs(p) + fabs(q) + fabs(r)) != 0.0) { | |
| p /= x; | |
| q /= x; | |
| r /= x; | |
| } | |
| } | |
| if ((s = SIGN(sqrt(p * p + q * q + r * r), p)) != 0.0) { | |
| if (k == m) { | |
| if (l != m) | |
| a[k][k - 1] = -a[k][k - 1]; | |
| } else | |
| a[k][k - 1] = -s * x; | |
| p += s; | |
| x = p / s; | |
| y = q / s; | |
| z = r / s; | |
| q /= p; | |
| r /= p; | |
| for (j = k; j < nn + 1; j++) { | |
| p = a[k][j] + q * a[k + 1][j]; | |
| if (k + 1 != nn) { | |
| p += r * a[k + 2][j]; | |
| a[k + 2][j] -= p * z; | |
| } | |
| a[k + 1][j] -= p * y; | |
| a[k][j] -= p * x; | |
| } | |
| mmin = nn < k + 3 ? nn : k + 3; | |
| for (i = l; i < mmin + 1; i++) { | |
| p = x * a[i][k] + y * a[i][k + 1]; | |
| if (k != (nn)) { | |
| p += z * a[i][k + 2]; | |
| a[i][k + 2] -= p * r; | |
| } | |
| a[i][k + 1] -= p * q; | |
| a[i][k] -= p; | |
| } | |
| } | |
| } | |
| } | |
| } | |
| } while (l + 1 < nn); | |
| } | |
| } | |
| /********************************************************* | |
| * calculate eigenvalues for a non-symmetric real matrix * | |
| *********************************************************/ | |
| void n_eigen(double *_a, int n, double *wr, double *wi) | |
| { | |
| int i; | |
| double **a = (double **) calloc(n, sizeof(void *)); | |
| for (i = 0; i < n; ++i) | |
| a[i] = _a + i * n; | |
| balanc(a, n); | |
| elmhes(a, n); | |
| hqr(a, n, wr, wi); | |
| free(a); | |
| } | |
| /* convert a symmetric matrix to tridiagonal form */ | |
| #define SQR(a) ((a)*(a)) | |
| static double pythag(double a, double b) | |
| { | |
| double absa, absb; | |
| absa = fabs(a); | |
| absb = fabs(b); | |
| if (absa > absb) return absa * sqrt(1.0 + SQR(absb / absa)); | |
| else return (absb == 0.0 ? 0.0 : absb * sqrt(1.0 + SQR(absa / absb))); | |
| } | |
| static void tred2(double **a, int n, double *d, double *e) | |
| { | |
| 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 < l + 1; k++) | |
| scale += fabs(a[i][k]); | |
| if (scale == 0.0) | |
| e[i] = a[i][l]; | |
| else { | |
| for (k = 0; k < l + 1; k++) { | |
| a[i][k] /= scale; | |
| h += a[i][k] * a[i][k]; | |
| } | |
| f = a[i][l]; | |
| g = (f >= 0.0 ? -sqrt(h) : sqrt(h)); | |
| e[i] = scale * g; | |
| h -= f * g; | |
| a[i][l] = f - g; | |
| f = 0.0; | |
| for (j = 0; j < l + 1; j++) { | |
| /* Next statement can be omitted if eigenvectors not wanted */ | |
| a[j][i] = a[i][j] / h; | |
| g = 0.0; | |
| for (k = 0; k < j + 1; k++) | |
| g += a[j][k] * a[i][k]; | |
| for (k = j + 1; k < l + 1; k++) | |
| g += a[k][j] * a[i][k]; | |
| e[j] = g / h; | |
| f += e[j] * a[i][j]; | |
| } | |
| hh = f / (h + h); | |
| for (j = 0; j < l + 1; j++) { | |
| f = a[i][j]; | |
| e[j] = g = e[j] - hh * f; | |
| for (k = 0; k < j + 1; k++) | |
| a[j][k] -= (f * e[k] + g * a[i][k]); | |
| } | |
| } | |
| } else | |
| e[i] = a[i][l]; | |
| d[i] = h; | |
| } | |
| /* Next statement can be omitted if eigenvectors not wanted */ | |
| d[0] = 0.0; | |
| e[0] = 0.0; | |
| /* Contents of this loop can be omitted if eigenvectors not wanted except for statement d[i]=a[i][i]; */ | |
| for (i = 0; i < n; i++) { | |
| l = i; | |
| if (d[i] != 0.0) { | |
| for (j = 0; j < l; j++) { | |
| g = 0.0; | |
| for (k = 0; k < l; k++) | |
| g += a[i][k] * a[k][j]; | |
| for (k = 0; k < l; k++) | |
| a[k][j] -= g * a[k][i]; | |
| } | |
| } | |
| d[i] = a[i][i]; | |
| a[i][i] = 1.0; | |
| for (j = 0; j < l; j++) | |
| a[j][i] = a[i][j] = 0.0; | |
| } | |
| } | |
| /* calculate the eigenvalues and eigenvectors of a symmetric tridiagonal matrix */ | |
| static void tqli(double *d, double *e, int n, double **z) | |
| { | |
| int m, l, iter, i, k; | |
| double s, r, p, g, f, dd, c, b; | |
| 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 = fabs(d[m]) + fabs(d[m + 1]); | |
| if (fabs(e[m]) + dd == dd) | |
| break; | |
| } | |
| if (m != l) { | |
| if (iter++ == 30) { | |
| fprintf(stderr, "[tqli] Too many iterations in tqli.\n"); | |
| break; | |
| } | |
| 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; | |
| /* Next loop can be omitted if eigenvectors not wanted */ | |
| 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); | |
| } | |
| } | |
| int n_eigen_symm(double *_a, int n, double *eval) | |
| { | |
| double **a, *e; | |
| int i; | |
| a = (double**)calloc(n, sizeof(void*)); | |
| e = (double*)calloc(n, sizeof(double)); | |
| for (i = 0; i < n; ++i) a[i] = _a + i * n; | |
| tred2(a, n, eval, e); | |
| tqli(eval, e, n, a); | |
| free(a); free(e); | |
| return 0; | |
| } | |
| #ifdef LH3_EIGEN_MAIN | |
| int main(void) | |
| { | |
| int i, j; | |
| static double u[5], v[5]; | |
| static double a[5][5] = {{1.0, 6.0, -3.0, -1.0, 7.0}, | |
| {8.0, -15.0, 18.0, 5.0, 4.0}, {-2.0, 11.0, 9.0, 15.0, 20.0}, | |
| {-13.0, 2.0, 21.0, 30.0, -6.0}, {17.0, 22.0, -5.0, 3.0, 6.0}}; | |
| static double b[5][5]={ {10.0,1.0,2.0,3.0,4.0}, | |
| {1.0,9.0,-1.0,2.0,-3.0}, | |
| {2.0,-1.0,7.0,3.0,-5.0}, | |
| {3.0,2.0,3.0,12.0,-1.0}, | |
| {4.0,-3.0,-5.0,-1.0,15.0}}; | |
| printf("MAT H IS:\n"); | |
| for (i = 0; i <= 4; i++) { | |
| for (j = 0; j <= 4; j++) | |
| printf("%13.7e ", a[i][j]); | |
| printf("\n"); | |
| } | |
| printf("\n"); | |
| n_eigen(a[0], 5, u, v); | |
| for (i = 0; i <= 4; i++) | |
| printf("%13.7e +J %13.7e\n", u[i], v[i]); | |
| printf("\n"); | |
| printf("\n\n"); | |
| n_eigen_symm(b[0], 5, u); | |
| for (i = 0; i <= 4; i++) | |
| printf("%13.7e\n", u[i]); | |
| printf("\n"); | |
| for (i = 0; i <= 4; i++) { | |
| for (j = 0; j <= 4; j++) | |
| printf("%12.6e ", b[i][j]); | |
| printf("\n"); | |
| } | |
| printf("\n"); | |
| return 0; | |
| } | |
| /* correct output: | |
| 4.2961008e+01 +J 0.0000000e+00 | |
| -6.6171383e-01 +J 0.0000000e+00 | |
| -1.5339638e+01 +J -6.7556929e+00 | |
| -1.5339638e+01 +J 6.7556929e+00 | |
| 1.9379983e+01 +J 0.0000000e+00 | |
| */ | |
| #endif |
Author
This was not written by me. I forgot where it comes from. I only remember that the original source code doesn't come with license information.
I appreciate your reply. For what it's worth, I finished testing my version and it's finally mature and safe to use now. (I finally added automated tests and benchmarks. Again, please forgive the self-promotion.) License issues are a huge headache.
Hello,
I was wondering if this process had a generic name. Like the power method can be used to find eigenvalue but this seems to be different for the symmetric matrix in particular. Any help on a source for the theory would be great!
p.s. code works great, i checked it with MATLABS eig() function and results are identical. (This function uses the power method fyi)
Also checkout this version: https://github.com/swedishembedded/control/blob/main/src/linalg/eig_sym.c
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Thanks
(License?)
(I have recently been having a hard time finding a simple eigenvector calculator with a permissive explicit license. Out of that frustration, I wrote this repository and released it in the public domain using the CC0 license. Please forgive the self-promotion.)