# include <float.h>
# include <math.h>
# include <stdio.h>
# include <stdlib.h>

# include "power_method.h"

/******************************************************************************/

void power_method ( int n, double a[], double y[], int it_max, double tol,
  double *lambda, int *it_num )

/******************************************************************************/
/*
  Purpose:

    power_method() applies several steps of the power method.

  Discussion:

    For a given NxN matrix A and an N vector Y, the power method produces
    a series of estimates for LAMBDA, the largest eigenvalue, and Y,
    the eigenvector corresponding to LAMBDA.

    The iteration repeats the following steps

      AY     = A * Y
      LAMBDA = || AY ||
      Y      = AY / LAMBDA

    If the matrix A has a single real eigenvalue of maximum modulus,
    then this iteration will generally produce a good estimate for that
    eigenvalue and its corresponding eigenvector.

    If there are multiple distinct eigenvalues of the same modulus,
    perhaps two values of opposite sign, or complex eigenvalues, then
    the situation is more complicated.

    Separate issues:

    * when estimating the value of LAMBDA, we use the Rayleigh quotient,
    LAMBDA = ( y' * A * y ) / ( y' * y ).  Since we normalize Y, the
    bottom of the fraction is 1.  Using this estimate allows us to
    easily capture the sign of LAMDBA.  Using the eucldean norm
    instead, for instance, would always give a positive value.

    * If the dominant eigenvalue is negative, then the iteration 
    as given will produce eigenvector iterates that alternate in sign.  
   
    * It is worth knowing whether the successive eigenvector estimates
    are tending to some value.  Since an eigenvector is really a direction,
    we need to normalize the vectors, and we need to somehow treat both
    a vector and its negative as holding the same information.  This
    means that the proper way of measuring the difference between two
    eigenvector estimates is to normalize them both, and then compute
    the cosine between them as y1'y2, followed by the sine, which is
    sqrt ( 1 - ( y1'y2)^2 ).  If this sine is small, the vectors y1 and y2
    are "close" in the sense of direction.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    18 July 2008

  Author:

    John Burkardt

  Input:

    int N, the order of the matrix.

    double A[N*N], the matrix.

    Input/double Y[N], the estimate for the eigenvector.

    int IT_MAX, the maximum number of iterations to take.
    1 <= IT_MAX.

    double TOL, an error tolerance.

  Output:

    double *LAMBDA, the estimate for the eigenvalue.

    int *IT_NUM, the number of iterations taken.
*/
{
  double *ay;
  double cos_y1y2;
  static int debug = 0;
  int i;
  int it;
  double lambda_old;
  double norm;
  double sin_y1y2;
  double val_dif;
  double *y_old;

  ay = ( double * ) malloc ( n * sizeof ( double ) );
  y_old = ( double * ) malloc ( n * sizeof ( double ) );

  if ( debug )
  {
    printf ( "\n" );
    printf ( "     IT      Lambda          Delta-Lambda    Delta-Y\n" );
    printf ( "\n" );
  }
/*
  Force Y to be a vector of unit norm.
*/
  norm = r8vec_norm_l2 ( n, y );

  for ( i = 0; i < n; i++ )
  {
    y[i] = y[i] / norm;
  }
  it = 0;

  for ( i = 0; i < n; i++ )
  {
    y_old[i] = y[i];
  }

  r8mat_mv ( n, n, a, y, ay );
  *lambda = r8vec_dot ( n, y, ay );
  norm = r8vec_norm_l2 ( n, ay );
  for ( i = 0; i < n; i++ )
  {
    y[i] = ay[i] / norm;
  }
  if ( *lambda < 0.0 )
  {
    for ( i = 0; i < n; i++ )
    {
      y[i] = - y[i];
    }
  }

  val_dif = 0.0;
  cos_y1y2 = r8vec_dot ( n, y, y_old );
  sin_y1y2 = sqrt ( ( 1.0 - cos_y1y2 ) * ( 1.0 + cos_y1y2 ) );

  if ( debug )
  {
    printf ( "  %5d  %14f  %14e  %14e\n", it, *lambda, val_dif, sin_y1y2 );
  }

  for ( it = 1; it <= it_max; it++ )
  {
    lambda_old = *lambda;
    for ( i = 0; i < n; i++ )
    {
      y_old[i] = y[i];
    }

    r8mat_mv ( n, n, a, y, ay );
    *lambda = r8vec_dot ( n, y, ay );
    norm = r8vec_norm_l2 ( n, ay );
    for ( i = 0; i < n; i++ )
    {
      y[i] = ay[i] / norm;
    }
    if ( *lambda < 0.0 )
    {
      for ( i = 0; i < n; i++ )
      {
        y[i] = - y[i];
      }
    }

    val_dif = fabs ( *lambda - lambda_old );
    cos_y1y2 = r8vec_dot ( n, y, y_old );
    sin_y1y2 = sqrt ( ( 1.0 - cos_y1y2 ) * ( 1.0 + cos_y1y2 ) );

    if ( debug )
    {
      printf ( "  %5d  %14f  %14e  %14e\n", it, *lambda, val_dif, sin_y1y2 );
    }

    if ( val_dif <= tol )
    {
      break;
    }
  }

  for ( i = 0; i < n; i++ )
  {
    y[i] = ay[i] / *lambda;
  }

  *it_num = it;

  free ( ay );
  free ( y_old );

  return;
}
/******************************************************************************/

void r8mat_mv ( int m, int n, double a[], double x[], double y[] )

/******************************************************************************/
/*
  Purpose:

    r8mat_mv() multiplies a matrix times a vector.

  Discussion: 							    

    An R8MAT is a doubly dimensioned array of double precision values, which
    may be stored as a vector in column-major order.

    For this routine, the result is returned as the function value.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 April 2007

  Author:

    John Burkardt

  Input:

    int M, N, the number of rows and columns of the matrix.

    double A[M,N], the M by N matrix.

    double X[N], the vector to be multiplied by A.

  Output:

    double Y[M], the product A*X.
*/
{
  int i;
  int j;

  for ( i = 0; i < m; i++ )
  {
    y[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      y[i] = y[i] + a[i+j*m] * x[j];
    }
  }

  return;
}
/******************************************************************************/

double r8vec_dot ( int n, double a1[], double a2[] )

/******************************************************************************/
/*
  Purpose:

    r8vec_dot() computes the dot product of a pair of R8VEC's.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    26 July 2007

  Author:

    John Burkardt

  Input:

    int N, the number of entries in the vectors.

    double A1[N], A2[N], the two vectors to be considered.

  Output:

    double R8VEC_DOT, the dot product of the vectors.
*/
{
  int i;
  double value;

  value = 0.0;
  for ( i = 0; i < n; i++ )
  {
    value = value + a1[i] * a2[i];
  }
  return value;
}
/******************************************************************************/

double r8vec_norm_l2 ( int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    r8vec_norm_l2() returns the L2 norm of an R8VEC.

  Discussion:

    The vector L2 norm is defined as:

      R8VEC_NORM_L2 = sqrt ( sum ( 1 <= I <= N ) A(I)**2 ).

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 March 2003

  Author:

    John Burkardt

  Input:

    int N, the number of entries in A.

    double A[N], the vector whose L2 norm is desired.

  Output:

    double R8VEC_NORM_L2, the L2 norm of A.
*/
{
  int i;
  double v;

  v = 0.0;

  for ( i = 0; i < n; i++ )
  {
    v = v + a[i] * a[i];
  }
  v = sqrt ( v );

  return v;
}