62c_vector<double, 3> CalculateEigenvectorForSmallestNonzeroEigenvalue(c_matrix<double, 3, 3>& rA)
65 if (norm_inf(rA - trans(rA)) > 10 * DBL_EPSILON)
73 c_matrix<double, 3, 3> copy_A(rA);
75 c_vector<double, 3> eigenvec1 = scalar_vector<double>(3, 1.0);
77 double eigen1 = CalculateMaxEigenpair(copy_A, eigenvec1);
80 c_matrix<double, 3, 3> wielandt_reduce_first_vector = identity_matrix<double>(3, 3);
81 wielandt_reduce_first_vector -= outer_prod(eigenvec1, eigenvec1);
82 copy_A = prod(wielandt_reduce_first_vector, copy_A);
84 c_vector<double, 3> eigenvec2 = scalar_vector<double>(3, 1.0);
85 double eigen2 = CalculateMaxEigenpair(copy_A, eigenvec2);
88 c_matrix<double, 3, 3> wielandt_reduce_second_vector = identity_matrix<double>(3, 3);
89 wielandt_reduce_second_vector -= outer_prod(eigenvec2, eigenvec2);
90 copy_A = prod(wielandt_reduce_second_vector, copy_A);
92 c_vector<double, 3> eigenvec3 = scalar_vector<double>(3, 1.0);
93 double eigen3 = CalculateMaxEigenpair(copy_A, eigenvec3);
96 if (eigen3 >= DBL_EPSILON)
100 if (eigen2 >= DBL_EPSILON)
105 assert(eigen1 > DBL_EPSILON);
109double CalculateMaxEigenpair(c_matrix<double, 3, 3>& rA, c_vector<double, 3>& rEigenvector)
112 double step = DBL_MAX;
113 while (step > DBL_EPSILON)
115 c_vector<double, 3> old_value(rEigenvector);
116 rEigenvector = prod(rA, rEigenvector);
117 norm = norm_2(rEigenvector);
118 rEigenvector /= norm;
119 if (norm < DBL_EPSILON)
124 step = norm_inf(rEigenvector - old_value);