Compadre  1.3.3
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 
6 namespace Compadre {
7 namespace GMLS_LinearAlgebra {
8 
9 KOKKOS_INLINE_FUNCTION
10 void 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] = {1,1,1};
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 
pool_type::generator_type generator_type
team_policy::member_type member_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.
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
Kokkos::Random_XorShift64_Pool pool_type