1 | |
---|
2 | #include <stdio.h> |
---|
3 | #include <math.h> |
---|
4 | |
---|
5 | void eigVl(float mat[3][3]) |
---|
6 | { |
---|
7 | |
---|
8 | float eigValue[3]; |
---|
9 | float eigVc[9]; |
---|
10 | |
---|
11 | float a = 0; |
---|
12 | float b = 0; |
---|
13 | |
---|
14 | float c[3]; |
---|
15 | |
---|
16 | // c[0] is the determinante of mat |
---|
17 | c[0] = mat[0][0] * mat[1][1] * mat[2][2] + |
---|
18 | 2* mat[0][1] * mat[0][2] * mat[1][2] - |
---|
19 | mat[0][0] * mat[1][2] * mat[1][2] - |
---|
20 | mat[1][1] * mat[0][2] * mat[0][2] - |
---|
21 | mat[2][2] * mat[0][1] * mat[0][1]; |
---|
22 | |
---|
23 | // c[1] is the trace of a |
---|
24 | c[1] = mat[0][0] * mat[1][1] - |
---|
25 | mat[0][1] * mat[0][1] + |
---|
26 | mat[0][0] * mat[2][2] - |
---|
27 | mat[0][2] * mat[0][2] + |
---|
28 | mat[1][1] * mat[2][2] - |
---|
29 | mat[1][2] * mat[1][2]; |
---|
30 | |
---|
31 | // c[2] is the sum of the diagonal elements |
---|
32 | c[2] = mat[0][0] + |
---|
33 | mat[1][1] + |
---|
34 | mat[2][2]; |
---|
35 | |
---|
36 | |
---|
37 | // Computing the roots: |
---|
38 | a = (3.0*c[1] - c[2]*c[2]) / 3.0; |
---|
39 | b = (-2.0*c[2]*c[2]*c[2] + 9.0*c[1]*c[2] - 27.0*c[0]) / 27.0; |
---|
40 | |
---|
41 | float Q = b*b/4.0 + a*a*a/27.0; |
---|
42 | |
---|
43 | if (Q < 0) |
---|
44 | { |
---|
45 | printf("good\n"); |
---|
46 | float psi = atan2(sqrt(-Q), -b/2.0); |
---|
47 | float p = sqrt((b/2.0)*(b/2.0) - Q); |
---|
48 | |
---|
49 | eigValue[0] = c[2]/3.0 + 2 * pow(p, 1/3.0) * cos(psi/3.0); |
---|
50 | eigValue[1] = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
51 | + sqrt(3.0) * sin(psi/3.0)); |
---|
52 | eigValue[2] = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
53 | - sqrt(3.0) * sin(psi/3.0)); |
---|
54 | |
---|
55 | } |
---|
56 | else if (Q == 0) |
---|
57 | { |
---|
58 | eigValue[0] = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
59 | eigValue[1] = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
60 | eigValue[2] = c[2]/3.0 + 2* pow(b/2.0, 1.0/3.0); |
---|
61 | } |
---|
62 | else if (Q > 0) |
---|
63 | { |
---|
64 | printf("A is multiple of Identity matrix (lambda * I3))\n"); |
---|
65 | eigValue[0] = eigValue[1] = eigValue[2] = 1; |
---|
66 | } |
---|
67 | |
---|
68 | printf("input: | %f | %f | %f |\n", mat[0][0], mat[0][1], mat[0][2] ); |
---|
69 | printf(" | %f | %f | %f |\n", mat[1][0], mat[1][1], mat[1][2] ); |
---|
70 | printf(" | %f | %f | %f |\n", mat[2][0], mat[2][1], mat[2][2] ); |
---|
71 | |
---|
72 | printf("%f %f %f\n", eigValue[0], eigValue[1], eigValue[2]); |
---|
73 | |
---|
74 | |
---|
75 | |
---|
76 | |
---|
77 | } |
---|