00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef EIGENVECTORS_HEADER
00030 #define EIGENVECTORS_HEADER
00031
00032
00033
00034
00035
00036
00047 #include "matrix_utilities.h"
00048 #include "integral_matrix_wrappers.h"
00049 #include "SizesAndBlocks.h"
00050 #include "output.h"
00051
00052 #include <iostream>
00053 #include <string.h>
00054
00055 #include "LanczosSeveralLargestEig.h"
00056
00057
00058 namespace eigvec
00059 {
00061 template<typename Treal, typename MatrixType, typename VectorType>
00062 Treal compute_rayleigh_quotient(const MatrixType& A, const VectorType& eigVec)
00063 {
00064 VectorType y, Ay;
00065
00066 y = eigVec;
00067 Treal ONE = (Treal)1.0;
00068 y *= ONE / y.eucl();
00069 Ay = ONE * A * y;
00070 Treal lambda = transpose(y) * Ay;
00071 return lambda;
00072 }
00073
00074
00077 template<typename Treal, typename MatrixType, typename VectorType>
00078 void lanczos_method(const MatrixType& A,
00079 std::vector<Treal>& eigVal,
00080 std::vector<VectorType>& eigVec,
00081 int number_of_eigenvalues,
00082 const Treal TOL,
00083 std::vector<int>& num_iter,
00084 int maxit = 200,
00085 bool do_deflation = false)
00086 {
00087 assert(eigVal.size() == eigVec.size());
00088 assert(eigVal.size() == num_iter.size());
00089 assert(number_of_eigenvalues >= 1);
00090
00091 if (eigVec[0].is_empty())
00092 {
00093 throw std::runtime_error("Error in eigvec::lanczos_method : eigVec[0].is_empty()");
00094 }
00095
00096 const Treal ONE = 1.0;
00097
00098 if (!do_deflation)
00099 {
00100 try
00101 {
00102 VectorType y;
00103 y = eigVec[0];
00104 y *= (ONE / y.eucl());
00105
00106
00107 mat::arn::LanczosSeveralLargestEig<Treal, MatrixType, VectorType> lan
00108 (A, y, number_of_eigenvalues, maxit);
00109 lan.setAbsTol(TOL);
00110 lan.setRelTol(TOL);
00111 lan.run();
00112 Treal acc = 0;
00113 lan.get_ith_eigenpair(1, eigVal[0], eigVec[0], acc);
00114
00115 VectorType resVec(eigVec[0]);
00116 resVec *= eigVal[0];
00117 resVec += -ONE * A * eigVec[0];
00118
00119
00120
00121
00122
00123
00124
00125
00126 num_iter[0] = lan.get_num_iter();
00127 }
00128 catch (std::exception& e)
00129 {
00130 num_iter[0] = maxit;
00131 }
00132 }
00133 else
00134 {
00135
00136 if (number_of_eigenvalues > 1)
00137 {
00138 VectorType y, v1;
00139 v1 = eigVec[0];
00140
00141
00142 if (eigVec[1].is_empty())
00143 {
00144 throw std::runtime_error("Error in eigvec::lanczos_method : eigVec[1].is_empty()");
00145 }
00146 y = eigVec[1];
00147 y *= (ONE / y.eucl());
00148
00149 try
00150 {
00151 int num_eig = 1;
00152
00153 Treal eigmin, eigmax;
00154 A.gershgorin(eigmin, eigmax);
00155 Treal sigma = eigVal[0] - eigmin;
00156 mat::arn::LanczosSeveralLargestEig<Treal, MatrixType, VectorType> lan
00157 (A, y, num_eig, maxit, 100, &v1, sigma);
00158 lan.setAbsTol(TOL);
00159 lan.setRelTol(TOL);
00160 lan.run();
00161 Treal acc = 0;
00162 lan.get_ith_eigenpair(1, eigVal[1], eigVec[1], acc);
00163
00164 VectorType resVec(eigVec[1]);
00165 resVec *= eigVal[1];
00166 resVec += -ONE * A * eigVec[1];
00167
00168 num_iter[1] = lan.get_num_iter();
00169 }
00170 catch (std::exception& e)
00171 {
00172 num_iter[1] = maxit;
00173 }
00174 }
00175 else
00176 {
00177 throw std::runtime_error("Error in eigvec::lanczos_method : number_of_eigenvalues <= 1");
00178 }
00179 }
00180 }
00181
00182
00185 template<typename Treal, typename MatrixType, typename VectorType>
00186 void power_method(const MatrixType& A,
00187 Treal& eigVal,
00188 VectorType& eigVec,
00189 const Treal TOL,
00190 int& num_iter,
00191 int maxit = 200)
00192 {
00193 VectorType y;
00194 VectorType Ay;
00195 VectorType residual;
00196 VectorType temp;
00197 Treal lambda;
00198 const Treal ONE = 1.0;
00199 const Treal MONE = -1.0;
00200
00201 y = eigVec;
00202 y *= (ONE / y.eucl());
00203
00204
00205 Ay = y;
00206 residual = y;
00207 temp = y;
00208
00209 int it = 1;
00210 Ay = ONE * A * y;
00211
00212 while (it == 1 || (residual.eucl() / template_blas_fabs(lambda) > TOL && it <= maxit))
00213 {
00214 y = Ay;
00215 y *= ONE / Ay.eucl();
00216 Ay = ONE * A * y;
00217 lambda = transpose(y) * Ay;
00218
00219
00220 residual = Ay;
00221 residual += (MONE * lambda) * y;
00222
00223
00224 ++it;
00225 }
00226
00227 printf("Power method required %d iterations.\n", it - 1);
00228
00229 eigVal = lambda;
00230 eigVec = y;
00231 num_iter = it - 1;
00232 }
00233
00234
00236 template<typename Treal, typename MatrixType, typename VectorType>
00237 int computeEigenvectors(const MatrixType& A,
00238 Treal tol,
00239 std::vector<Treal>& eigVal,
00240 std::vector<VectorType>& eigVec,
00241 int number_of_eigenvalues_to_compute,
00242 std::string method,
00243 std::vector<int>& num_iter,
00244 int maxit = 200,
00245 bool do_deflation = false
00246 )
00247 {
00248 assert(number_of_eigenvalues_to_compute >= 1);
00249 assert(eigVal.size() >= 1);
00250 assert(eigVec.size() == eigVal.size());
00251 assert(eigVec.size() == num_iter.size());
00252
00253 if (method == "power")
00254 {
00255 if (eigVal.size() > 1)
00256 {
00257 throw std::runtime_error("Error in eigvec::computeEigenvectors: computation of more "
00258 "than 1 eigenpair is not implemented for the power method.");
00259 }
00260 if (do_deflation)
00261 {
00262 throw std::runtime_error("Error in eigvec::computeEigenvectors: deflation is not implemented for the power method.");
00263 }
00264 power_method(A, eigVal[0], eigVec[0], tol, num_iter[0], maxit);
00265 }
00266 else if (method == "lanczos")
00267 {
00268 lanczos_method(A, eigVal, eigVec, number_of_eigenvalues_to_compute, tol, num_iter, maxit, do_deflation);
00269 }
00270 else
00271 {
00272 throw std::runtime_error("Error in eigvec::computeEigenvectors: unknown method.");
00273 }
00274 return 0;
00275 }
00276 }
00277
00278 #endif // EIGENVECTORS_HEADER