29 #ifndef EIGENVECTORS_HEADER
30 #define EIGENVECTORS_HEADER
61 template<
typename Treal,
typename MatrixType,
typename VectorType>
67 Treal ONE = (Treal)1.0;
77 template<
typename Treal,
typename MatrixType,
typename VectorType>
79 std::vector<Treal>& eigVal,
80 std::vector<VectorType>& eigVec,
81 int number_of_eigenvalues,
83 std::vector<int>& num_iter,
85 bool do_deflation =
false)
87 assert(eigVal.size() == eigVec.size());
88 assert(eigVal.size() == num_iter.size());
89 assert((
int) eigVal.size() >= number_of_eigenvalues);
90 assert(number_of_eigenvalues >= 1);
92 if (eigVec[0].is_empty())
94 throw std::runtime_error(
"Error in eigvec::lanczos_method : eigVec[0].is_empty()");
97 const Treal ONE = 1.0;
103 y *= (ONE / y.
eucl());
106 (
A, y, number_of_eigenvalues, maxit);
121 for (
int i = 1; i <= number_of_eigenvalues; i++) {
130 if (number_of_eigenvalues > 1)
136 if (eigVec[1].is_empty())
138 throw std::runtime_error(
"Error in eigvec::lanczos_method : eigVec[1].is_empty()");
141 y *= (ONE / y.
eucl());
147 Treal eigmin, eigmax;
148 A.gershgorin(eigmin, eigmax);
149 Treal sigma = eigVal[0] - eigmin;
151 (
A, y, num_eig, maxit, 100, &v1, sigma);
160 resVec += -ONE *
A * eigVec[1];
172 throw std::runtime_error(
"Error in eigvec::lanczos_method : number_of_eigenvalues <= 1");
180 template<
typename Treal,
typename MatrixType,
typename VectorType>
193 const Treal ONE = 1.0;
194 const Treal MONE = -1.0;
197 y *= (ONE / y.
eucl());
210 y *= ONE / Ay.
eucl();
216 residual += (MONE * lambda) * y;
222 printf(
"Power method required %d iterations.\n", it - 1);
231 template<
typename Treal,
typename MatrixType,
typename VectorType>
234 std::vector<Treal>& eigVal,
235 std::vector<VectorType>& eigVec,
236 int number_of_eigenvalues_to_compute,
238 std::vector<int>& num_iter,
240 bool do_deflation =
false
243 assert(number_of_eigenvalues_to_compute >= 1);
244 assert(eigVal.size() >= 1);
245 assert(eigVec.size() == eigVal.size());
246 assert(eigVec.size() == num_iter.size());
248 if (method ==
"power")
250 if (eigVal.size() > 1)
252 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: computation of more "
253 "than 1 eigenpair is not implemented for the power method.");
257 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: deflation is not implemented for the power method.");
259 power_method(
A, eigVal[0], eigVec[0], tol, num_iter[0], maxit);
261 else if (method ==
"lanczos")
263 lanczos_method(
A, eigVal, eigVec, number_of_eigenvalues_to_compute, tol, num_iter, maxit, do_deflation);
267 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: unknown method.");
273 #endif // EIGENVECTORS_HEADER