Compadre 1.5.5
Loading...
Searching...
No Matches
Compadre_LinearAlgebra_Definitions.hpp
Go to the documentation of this file.
1#ifndef _COMPADRE_LINEAR_ALGEBRA_DEFINITIONS_HPP_
2#define _COMPADRE_LINEAR_ALGEBRA_DEFINITIONS_HPP_
3
5
6namespace Compadre {
7namespace GMLS_LinearAlgebra {
8
9KOKKOS_INLINE_FUNCTION
10void largestTwoEigenvectorsThreeByThreeSymmetric(const member_type& teamMember, scratch_matrix_right_type V, scratch_matrix_right_type PtP, const int dimensions, pool_type& random_number_pool) {
11
12 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
13
14 double maxRange = 100;
15
16 generator_type rand_gen = random_number_pool.get_state();
17 // put in a power method here and a deflation by first found eigenvalue
18 double eigenvalue_relative_tolerance = 1e-6; // TODO: use something smaller, but really anything close is acceptable for this manifold
19
20
21 double v[3] = {rand_gen.drand(maxRange),rand_gen.drand(maxRange),rand_gen.drand(maxRange)};
22 double v_old[3] = {v[0], v[1], v[2]};
23
24 double error = 1;
25 double norm_v;
26
27 while (error > eigenvalue_relative_tolerance) {
28
29 double tmp1 = v[0];
30 v[0] = PtP(0,0)*tmp1 + PtP(0,1)*v[1];
31 if (dimensions>2) v[0] += PtP(0,2)*v[2];
32
33 double tmp2 = v[1];
34 v[1] = PtP(1,0)*tmp1 + PtP(1,1)*tmp2;
35 if (dimensions>2) v[1] += PtP(1,2)*v[2];
36
37 if (dimensions>2)
38 v[2] = PtP(2,0)*tmp1 + PtP(2,1)*tmp2 + PtP(2,2)*v[2];
39
40 norm_v = v[0]*v[0] + v[1]*v[1];
41 if (dimensions>2) norm_v += v[2]*v[2];
42 norm_v = std::sqrt(norm_v);
43
44 v[0] = v[0] / norm_v;
45 v[1] = v[1] / norm_v;
46 if (dimensions>2) v[2] = v[2] / norm_v;
47
48 error = (v[0]-v_old[0])*(v[0]-v_old[0]) + (v[1]-v_old[1])*(v[1]-v_old[1]);
49 if (dimensions>2) error += (v[2]-v_old[2])*(v[2]-v_old[2]);
50 error = std::sqrt(error);
51 error /= norm_v;
52
53
54 v_old[0] = v[0];
55 v_old[1] = v[1];
56 if (dimensions>2) v_old[2] = v[2];
57 }
58
59 double dot_product;
60 double norm;
61
62 // if 2D, orthonormalize second vector
63 if (dimensions==2) {
64
65 for (int i=0; i<2; ++i) {
66 V(0,i) = v[i];
67 }
68
69 // orthonormalize second eigenvector against first
70 V(1,0) = 1.0; V(1,1) = 1.0;
71 dot_product = V(0,0)*V(1,0) + V(0,1)*V(1,1);
72 V(1,0) -= dot_product*V(0,0);
73 V(1,1) -= dot_product*V(0,1);
74
75 norm = std::sqrt(V(1,0)*V(1,0) + V(1,1)*V(1,1));
76 V(1,0) /= norm;
77 V(1,1) /= norm;
78
79 } else { // otherwise, work on second eigenvalue
80
81 for (int i=0; i<3; ++i) {
82 V(0,i) = v[i];
83 for (int j=0; j<3; ++j) {
84 PtP(i,j) -= norm_v*v[i]*v[j];
85 }
86 }
87
88 error = 1;
89 v[0] = rand_gen.drand(maxRange); v[1] = rand_gen.drand(maxRange); v[2] = rand_gen.drand(maxRange);
90 v_old[0] = v[0]; v_old[1] = v[1]; v_old[2] =v[2];
91 while (error > eigenvalue_relative_tolerance) {
92
93 double tmp1 = v[0];
94 v[0] = PtP(0,0)*tmp1 + PtP(0,1)*v[1] + PtP(0,2)*v[2];
95
96 double tmp2 = v[1];
97 v[1] = PtP(1,0)*tmp1 + PtP(1,1)*tmp2 + PtP(1,2)*v[2];
98
99 v[2] = PtP(2,0)*tmp1 + PtP(2,1)*tmp2 + PtP(2,2)*v[2];
100
101 norm_v = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
102
103 v[0] = v[0] / norm_v;
104 v[1] = v[1] / norm_v;
105 v[2] = v[2] / norm_v;
106
107 error = std::sqrt((v[0]-v_old[0])*(v[0]-v_old[0]) + (v[1]-v_old[1])*(v[1]-v_old[1]) + (v[2]-v_old[2])*(v[2]-v_old[2])) / norm_v;
108
109 v_old[0] = v[0];
110 v_old[1] = v[1];
111 v_old[2] = v[2];
112 }
113
114 for (int i=0; i<3; ++i) {
115 V(1,i) = v[i];
116 }
117
118 // orthonormalize second eigenvector against first
119 dot_product = V(0,0)*V(1,0) + V(0,1)*V(1,1) + V(0,2)*V(1,2);
120
121 V(1,0) -= dot_product*V(0,0);
122 V(1,1) -= dot_product*V(0,1);
123 V(1,2) -= dot_product*V(0,2);
124
125 norm = std::sqrt(V(1,0)*V(1,0) + V(1,1)*V(1,1) + V(1,2)*V(1,2));
126 V(1,0) /= norm;
127 V(1,1) /= norm;
128 V(1,2) /= norm;
129
130 // get normal from cross product
131 V(2,0) = V(0,1)*V(1,2) - V(1,1)*V(0,2);
132 V(2,1) = V(1,0)*V(0,2) - V(0,0)*V(1,2);
133 V(2,2) = V(0,0)*V(1,1) - V(1,0)*V(0,1);
134
135 // orthonormalize third eigenvector (just to be sure)
136 norm = std::sqrt(V(2,0)*V(2,0) + V(2,1)*V(2,1) + V(2,2)*V(2,2));
137 V(2,0) /= norm;
138 V(2,1) /= norm;
139 V(2,2) /= norm;
140
141 }
142
143 random_number_pool.free_state(rand_gen);
144 });
145
146}
147
148} // GMLS_LinearAlgebra
149} // Compadre
150
151#endif
152
Kokkos::Random_XorShift64_Pool pool_type
team_policy::member_type member_type
pool_type::generator_type generator_type
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
KOKKOS_INLINE_FUNCTION void largestTwoEigenvectorsThreeByThreeSymmetric(const member_type &teamMember, scratch_matrix_right_type V, scratch_matrix_right_type PtP, const int dimensions, pool_type &random_number_pool)
Calculates two eigenvectors corresponding to two dominant eigenvalues.