1 #ifndef _COMPADRE_LINEAR_ALGEBRA_DEFINITIONS_HPP_ 2 #define _COMPADRE_LINEAR_ALGEBRA_DEFINITIONS_HPP_ 7 namespace GMLS_LinearAlgebra {
12 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
14 double maxRange = 100;
18 double eigenvalue_relative_tolerance = 1e-6;
21 double v[3] = {1,1,1};
22 double v_old[3] = {v[0], v[1], v[2]};
27 while (error > eigenvalue_relative_tolerance) {
30 v[0] = PtP(0,0)*tmp1 + PtP(0,1)*v[1];
31 if (dimensions>2) v[0] += PtP(0,2)*v[2];
34 v[1] = PtP(1,0)*tmp1 + PtP(1,1)*tmp2;
35 if (dimensions>2) v[1] += PtP(1,2)*v[2];
38 v[2] = PtP(2,0)*tmp1 + PtP(2,1)*tmp2 + PtP(2,2)*v[2];
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);
46 if (dimensions>2) v[2] = v[2] / norm_v;
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);
56 if (dimensions>2) v_old[2] = v[2];
65 for (
int i=0; i<2; ++i) {
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);
75 norm = std::sqrt(V(1,0)*V(1,0) + V(1,1)*V(1,1));
81 for (
int i=0; i<3; ++i) {
83 for (
int j=0; j<3; ++j) {
84 PtP(i,j) -= norm_v*v[i]*v[j];
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) {
94 v[0] = PtP(0,0)*tmp1 + PtP(0,1)*v[1] + PtP(0,2)*v[2];
97 v[1] = PtP(1,0)*tmp1 + PtP(1,1)*tmp2 + PtP(1,2)*v[2];
99 v[2] = PtP(2,0)*tmp1 + PtP(2,1)*tmp2 + PtP(2,2)*v[2];
101 norm_v = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
103 v[0] = v[0] / norm_v;
104 v[1] = v[1] / norm_v;
105 v[2] = v[2] / norm_v;
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;
114 for (
int i=0; i<3; ++i) {
119 dot_product = V(0,0)*V(1,0) + V(0,1)*V(1,1) + V(0,2)*V(1,2);
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);
125 norm = std::sqrt(V(1,0)*V(1,0) + V(1,1)*V(1,1) + V(1,2)*V(1,2));
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);
136 norm = std::sqrt(V(2,0)*V(2,0) + V(2,1)*V(2,1) + V(2,2)*V(2,2));
143 random_number_pool.free_state(rand_gen);
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