| 46 | } |
| 47 | |
| 48 | inline void Tred2(Mat3d &matrix_v, Vec3d &vector_d, Vec3d &vector_e) { |
| 49 | |
| 50 | // This is derived from the Algol procedures tred2 by |
| 51 | // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for |
| 52 | // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding |
| 53 | // Fortran subroutine in EISPACK. |
| 54 | |
| 55 | int i, j, k; |
| 56 | double f, g, h, hh; |
| 57 | for (j = 0; j < 3; j++) { |
| 58 | vector_d(j) = matrix_v(2, j); |
| 59 | } |
| 60 | |
| 61 | // Householder reduction to tridiagonal form. |
| 62 | |
| 63 | for (i = 3 - 1; i > 0; i--) { |
| 64 | |
| 65 | // Scale to avoid under/overflow. |
| 66 | |
| 67 | double scale = 0.0; |
| 68 | double h = 0.0; |
| 69 | for (k = 0; k < i; k++) { |
| 70 | scale = scale + fabs(vector_d(k)); |
| 71 | } |
| 72 | if (scale == 0.0) { |
| 73 | vector_e[i] = vector_d(i - 1); |
| 74 | for (j = 0; j < i; j++) { |
| 75 | vector_d(j) = matrix_v(i - 1, j); |
| 76 | matrix_v(i, j) = 0.0; |
| 77 | matrix_v(j, i) = 0.0; |
| 78 | } |
| 79 | } else { |
| 80 | |
| 81 | // Generate Householder vector. |
| 82 | |
| 83 | for (k = 0; k < i; k++) { |
| 84 | vector_d(k) /= scale; |
| 85 | h += vector_d(k) * vector_d(k); |
| 86 | } |
| 87 | f = vector_d(i - 1); |
| 88 | g = sqrt(h); |
| 89 | if (f > 0) { |
| 90 | g = -g; |
| 91 | } |
| 92 | vector_e(i) = scale * g; |
| 93 | h = h - f * g; |
| 94 | vector_d(i - 1) = f - g; |
| 95 | for (j = 0; j < i; j++) { |
| 96 | vector_e(j) = 0.0; |
| 97 | } |
| 98 | |
| 99 | // Apply similarity transformation to remaining columns. |
| 100 | |
| 101 | for (j = 0; j < i; j++) { |
| 102 | f = vector_d(j); |
| 103 | matrix_v(j, i) = f; |
| 104 | g = vector_e(j) + matrix_v(j, j) * f; |
| 105 | for (k = j + 1; k <= i - 1; k++) { |
no outgoing calls
no test coverage detected