14 #include <cuda/CudaGlobal.hpp>
23 #define MAX(a, b) ((a) > (b) ? (a) : (b))
28 static double hypot2(
double x,
double y) {
29 return sqrt(x * x + y * y);
35 static void tred2(
double V[n][n],
double d[n],
double e[n])
43 for(
int j = 0; j < n; j++) {
49 for(
int i = n - 1; i > 0; i--) {
55 for(
int k = 0; k < i; k++) {
56 scale = scale +
fabs(d[k]);
60 for(
int j = 0; j < i; j++) {
69 for(
int k = 0; k < i; k++) {
81 for(
int j = 0; j < i; j++) {
87 for(
int j = 0; j < i; j++) {
90 g = e[j] + V[j][j] * f;
91 for(
int k = j + 1; k <= i - 1; k++) {
98 for(
int j = 0; j < i; j++) {
102 double hh = f / (h + h);
103 for(
int j = 0; j < i; j++) {
106 for(
int j = 0; j < i; j++) {
109 for(
int k = j; k <= i - 1; k++) {
110 V[k][j] -= (f * e[k] + g * d[k]);
121 for(
int i = 0; i < n - 1; i++) {
122 V[n - 1][i] = V[i][i];
126 for(
int k = 0; k <= i; k++) {
127 d[k] = V[k][i + 1] / h;
129 for(
int j = 0; j <= i; j++) {
131 for(
int k = 0; k <= i; k++) {
132 g += V[k][i + 1] * V[k][j];
134 for(
int k = 0; k <= i; k++) {
139 for(
int k = 0; k <= i; k++) {
143 for(
int j = 0; j < n; j++) {
147 V[n - 1][n - 1] = 1.0;
153 static void tql2(
double V[n][n],
double d[n],
double e[n])
161 for(
int i = 1; i < n; i++) {
168 double eps = pow(2.0, -52.0);
169 for(
int l = 0; l < n; l++) {
173 tst1 = MAX(tst1,
fabs(d[l]) +
fabs(e[l]));
176 if(
fabs(e[m]) <= eps * tst1) {
193 double p = (d[l + 1] - g) / (2.0 * e[l]);
194 double r = hypot2(p, 1.0);
198 d[l] = e[l] / (p + r);
199 d[l + 1] = e[l] * (p + r);
200 double dl1 = d[l + 1];
202 for(
int i = l + 2; i < n; i++) {
213 double el1 = e[l + 1];
216 for(
int i = m - 1; i >= l; i--) {
226 p = c * d[i] - s * g;
227 d[i + 1] = h + s * (c * g + s * d[i]);
231 for(
int k = 0; k < n; k++) {
233 V[k][i + 1] = s * V[k][i] + c * h;
234 V[k][i] = c * V[k][i] - s * h;
237 p = -s * s2 * c3 * el1 * e[l] / dl1;
243 }
while(
fabs(e[l]) > eps * tst1);
251 for(
int i = 0; i < n - 1; i++) {
254 for(
int j = i + 1; j < n; j++) {
263 for(
int j = 0; j < n; j++) {
274 void EigenDecompSym3x3(
double A[n][n],
double V[n][n],
double d[n])
277 for(
int i = 0; i < n; i++) {
278 for(
int j = 0; j < n; j++) {
CU_HOST_DEVICE Vector< dim, T > fabs(const Vector< dim, T > &v)
Return the vector whose component is the absolute value of the input vector.
Definition: Vector.hpp:1175