/** @(#) numericalap/eigenpower.cpp */ #include "eigenpower.h" /** Realiza M pasos iterativos del metodo de la potencia * para obtener el eigenvalor dominante de la matriz A. * Tambien incorpora un proceso de ortogonalizacion * para encontrar los sucesivos eigenvalores. * @author Daniel Alba Cuellar, Omar Posada Villarreal (generalizar) * @param M numero de iteraciones * @param seq se intenta encontrar el eigenvalor numero seq * @return Maximo eigen valor. */ double eigenPower(SquaredMatrix &A, int M, int seq) { int k, i, j, l; double r; // eigen valor, lambda double max_y; VectorDouble y[MAX_DIM]; VectorDouble v_temp[MAX_DIM]; /* vector temporal */ double e_temp; /* escalar temporal */ /* w = x */ for (j=1; j <= n; j++) w[j] = x[j]; /* Proceso de ortogonalizacion: construye un vector inicial w a ser usado por el metodo de la potencia */ /* w = x - suma((V(.,j)*x/V(.,j)*V(.,j))V(.,j)). V(.,i) es el j-ésimo eigenvector de A. */ for (j=1; j < seq; j++) { /* extrae el j-ésimo eigenvector de la matriz V */ for (l=1; l <= n; l++) v_temp[l] = V[l][j]; e_temp = VtPorW(v_temp,x)/VtPorW(v_temp,v_temp); for (l=1; l <= n; l++) w[l] = w[l] - e_temp * v_temp[l]; } /* aqui comienza el metodo de la potencia */ for (k = 1; k <= M; k++) { /* y = Aw */ MporV( A, w, y ); /* r = w*y/w*w. esta es la aproximacion del eigenvalor por medio del cociente de Rayleigh */ r = VtPorW( w, y )/VtPorW( w, w ); /* calcula la norma infinito de y */ max_y = fabs(y[1]); for (i = 2; i <= n; i++) { if (fabs(y[i]) > max_y) { max_y = fabs(y[i]); } } /* w = y/norma(y). norma(y) es la norma infinito de y */ for (i = 1; i <= n; i++) { y[i] = y[i]/max_y; w[i] = y[i]; } /* Proceso de ortogonalizacion: elimina de w los errores de redondeo */ /* w = x - suma((V(.,j)*x/V(.,j)*V(.,j))V(.,j)). V(.,i) es el j-ésimo eigenvector de A. */ for (j=1; j < seq; j++) { /* extrae el j-ésimo eigenvector de la matriz V */ for (l=1; l <= n; l++) v_temp[l] = V[l][j]; e_temp = VtPorW(v_temp,y)/VtPorW(v_temp,v_temp); for (l=1; l <= n; l++) w[l] = w[l] - e_temp * v_temp[l]; } } imprimeW(); return r; } // Fin-------------------------------------------------------------------------