/** @(#) 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-------------------------------------------------------------------------