/* eigen_power.c   works sometimes, slow              */
/*          just finds the largest eigenvalue         */
/*                 given A, guess V = all ones        */
/*          loop:  V2 = A * V                         */
/*                 V3  = V2/largest_magnitude_element */
/*                 stop when |V3-V| small             */
/*                 V = V3   goto loop                 */
/*                                                    */
/*          check the eigenvalue and eigenvector      */

#include <stdio.h>
#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
static double eigenp(int n, double A[n][n], double V[], double *e);
static vecprint(int n, double V[]);
static matprint(int n, double A[][]);
static matmul(int n, double A[n][n], double V[], double V2[]);
static double determinant(int n, double A[n][n]);
static void inverse(int n, double A[n][n], double AA[n][n]);

int main(int argc, char *argv[])
{
  int i, j, k, n;
  double e, err, det;

  double A1[3][3] = {{ 3.0,  -1.0,   0.0},
                     {-2.0,   4.0,  -3.0},
                     { 0.0,  -1.0,   1.0}};

  double A2[3][3] = {{ 6.0,  12.0,  19.0},
                     {-9.0, -20.0, -33.0},
                     {  4.0,  9.0,  15.0}};

  double V[3], Vc[3], A[3][3], A3[3][3];

  printf("eigen_power.c running, find largest eigenvalue \n");
  n = 3;

  printf("initial matrix A1 n=%d\n", n);
  matprint(n, A1);

  err = eigenp(n, A1, V, &e);
  printf("eigenvalue=%12.8f, error=%g, eigenvector is \n", e, err);
  vecprint(n, V);

  printf("check eigenvalue and eigenvector \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      A[i][j] = A1[i][j];
      if(i==j) A[i][j] -= e;
    }
  }
  det = determinant(n, A);
  printf("determinant should be zero = %g \n", det);
  /* check A1 V=e V */
  for(i=0; i<n; i++)
  {
    Vc[i] = -e*V[i]; /* should zero */
    for(j=0; j<n; j++)
    {
      Vc[i]=Vc[i]+A1[i][j]*V[j];
    }
  }
  printf("eigenvector errors \n");
  vecprint(n, Vc);


  printf("\ninitial matrix A2 n=%d\n", n);
  matprint(n, A2);

  err = eigenp(n, A2, V, &e);
  printf("eigenvalue=%12.8f, error=%g, eigenvector is \n", e, err);
  vecprint(n, V);

  printf("check eigenvalue and eigenvector \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      A[i][j] = A2[i][j];
      if(i==j) A[i][j] -= e;
    }
  }
  det = determinant(n, A);
  printf("determinant should be zero = %g \n", det);
  /* check A2 V=e V */
  for(i=0; i<n; i++)
  {
    Vc[i] = -e*V[i]; /* should zero */
    for(j=0; j<n; j++)
    {
      Vc[i]=Vc[i]+A2[i][j]*V[j];
    }
  }
  printf("eigenvector errors \n");
  vecprint(n, Vc);

  printf("\ninvert A1 and find the smallest eigenvalue \n");
  inverse(n, A1, A3);
  printf("initial matrix A3 n=%d\n", n);
  matprint(n, A3);

  err = eigenp(n, A3, V, &e);
  printf("eigenvalue=%12.8f, error=%g, eigenvector is \n", e, err);
  vecprint(n, V);

  printf("check eigenvalue and eigenvector \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      A[i][j] = A3[i][j];
      if(i==j) A[i][j] -= e;
    }
  }
  det = determinant(n, A);
  printf("determinant should be zero = %g \n", det);
  /* check A3 V=e V */
  for(i=0; i<n; i++)
  {
    Vc[i] = -e*V[i]; /* should zero */
    for(j=0; j<n; j++)
    {
      Vc[i]=Vc[i]+A3[i][j]*V[j];
    }
  }
  printf("eigenvector errors \n");
  vecprint(n, Vc);
  printf("The smallest eigenvalue of A1 is 1 over, = %12.8f \n\n",
         1.0/e);
  printf("check eigenvalue of original A1 \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      A[i][j] = A1[i][j];
      if(i==j) A[i][j] -= 1.0/e;
    }
  }
  det = determinant(n, A);
  printf("determinant should be zero = %g \n", det);


  return 0;
} /* end main of eigen_power.c */

static double eigenp(int n, double A[n][n], double V[], double *e)
{
  int i, k;
  double big, bigabs, err;
  double V2[100], V3[100];
  int debug = 0;

  for(i=0; i<n; i++) V[i]=1.0;
  if(debug) {printf("initial guess V \n"); vecprint(n, V);}
  for(k=1; k<200; k++)
  {
    matmul(n, A, V, V2);
    if(debug) {printf("V2 = A * V \n"); vecprint(n, V2);}
    big = V2[0];
    bigabs = abs(big);
    for(i=1; i<n; i++)
    {
      if(abs(V2[i])>bigabs) big=V2[i];
      bigabs = abs(big);
    }
    for(i=0; i<n; i++) V3[i] = V2[i]/big;
    if(debug) {printf("normalized V3 \n"); vecprint(n, V3);}
    err = 0.0;
    for(i=0; i<n; i++) err = err + abs(V[i]- V3[i]);
    if(debug) printf("iteration k=%d, err=%g \n", k, err);
    for(i=0; i<n; i++) V[i] = V3[i];
    if(err<1.0e-6) break;
  }
  *e = big;
  return err;
} /* end eigenp */

static vecprint(int n, double V[])
{
  int i;
  for(i=0; i<n; i++)
  {
    printf("%12.8f  ", V[i]);
  }
  printf("\n\n");
} /* end vecprint */

static matprint(int n, double A[n][n])
{
  int i, j;
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      printf("%12.8f  ", A[i][j]);
    }
    printf("\n");
  }
  printf("\n");
} /* end matprint */

static matmul(int n, double A[n][n], double V[], double V2[])
{
  int i, j;
  for(i=0; i<n; i++)
  {
    V2[i] = 0.0;
    for(j=0; j<n; j++)
    {
       V2[i] = V2[i] + A[i][j]*V[j];
    }
  }
} /* end matmul */

static double determinant(int n, double A[n][n])
{
    double D = 1.0;  /* DETERMINANT */
    double *B;  /* WORKING MATRIX */
    int *ROW;  /* ROW INTERCHANGE INDICIES */
    int HOLD , I_PIVOT;  /* PIVOT INDICIES */
    double PIVOT;  /* PIVOT ELEMENT VALUE */
    double ABS_PIVOT;
    int i, j, k;

    B = (double *)calloc(n*n, sizeof(double));
    ROW = (int *)calloc(n, sizeof(int));
    memcpy(B,A,n*n*sizeof(double));

    /* SET UP ROW  INTERCHANGE VECTORS */
    for(k=0; k<n; k++){
      ROW[k]= k;
    }

    /* BEGIN MAIN REDUCTION LOOP  */
    for(k=0; k<n-1; k++){

      /* FIND LARGEST ELEMENT FOR PIVOT */
      PIVOT = B[ROW[k]*n+k];
      ABS_PIVOT = abs ( PIVOT );
      I_PIVOT = k;
      for(i=k; i<n; i++){
        if( abs(B[ROW[i]*n+k]) > ABS_PIVOT ){
          I_PIVOT = i;
          PIVOT = B[ROW[i]*n+k];
          ABS_PIVOT = abs(PIVOT);
        }
      }

      /* HAVE PIVOT, INTERCHANGE ROW POINTERS */
      if(I_PIVOT != k){
        HOLD = ROW[k];
        ROW[k] = ROW[I_PIVOT];
        ROW[I_PIVOT] = HOLD;
        D = - D;
      }

      /* CHECK FOR NEAR SINGULAR  */
      if(ABS_PIVOT < 1.0E-10){
        free(B);
        free(ROW);
        return 0.0;
      }
      else
      {
        D = D * PIVOT;

        /* REDUCE ABOUT PIVOT */
        for(j=k+1; j<n; j++){
          B[ROW[k]*n+j] = B[ROW[k]*n+j] / B[ROW[k]*n+k];
        }

        /* INNER REDUCTION LOOP */
        for(i=0; i<n; i++){
          if(i != k){
            for(j=k+1; j<n; j++){
              B[ROW[i]*n+j] = B[ROW[i]*n+j] - B[ROW[i]*n+k]* B[ROW[k]*n+j];
            }
          }
        }
      }
      /* FINISHED INNER REDUCTION */
    }
    /* END OF MAIN REDUCTION LOOP */
    return D * B[ROW[n-1]*n+n-1];
} /* end determinant */

static void inverse(int n, double A[n][n], double AA[n][n])
{
  int *row, *col;
  double *temp;
  int hold , I_pivot , J_pivot;
  double pivot, abs_pivot;
  int i, j, k;

  row = (int *)calloc(n, sizeof(int));
  col = (int *)calloc(n, sizeof(int));
  temp = (double *)calloc(n, sizeof(double));
  memcpy(AA,A,n*n*sizeof(double));

  /* set up row and column interchange vectors */
  for (k=0; k<n; k++)
  {
    row[k] = k ;
    col[k] = k ;
  }

  /* begin main reduction loop */
  for (k=0; k<n; k++) {

    /* find largest element for pivot */
    pivot = AA[row[k]][col[k]] ;
    I_pivot = k;
    J_pivot = k;
    for(i=k; i<n; i++){
      for(j=k; j<n; j++){
        abs_pivot = abs(pivot) ;
        if( abs(AA[row[i]][col[j]]) > abs_pivot ){
          I_pivot = i ;
          J_pivot = j ;
          pivot = AA[row[i]][col[j]] ;
        }
      }
    }
    if(abs(pivot) < 1.0E-10){
      free(row);
      free(col);
      free(temp);
      printf("MATRIX is singular !!! \n");
      return;
    }

    hold = row[k];
    row[k]= row[I_pivot];
    row[I_pivot] = hold ;
    hold = col[k];
    col[k]= col[J_pivot];
    col[J_pivot] = hold ;

    /* reduce about pivot */
    AA[row[k]][col[k]] = 1.0 / pivot ;
    for (j=0; j<n; j++) {
      if( j != k ){
        AA[row[k]][col[j]] = AA[row[k]][col[j]] * AA[row[k]][col[k]];
      }
    }

    /* inner reduction loop */
    for (i=0; i<n; i++)
    {
      if( k != i )
      {
        for (j=0; j<n; j++)
        {
          if( k != j )
          {
            AA[row[i]][col[j]] = AA[row[i]][col[j]] - AA[row[i]][col[k]] *
                                 AA[row[k]][col[j]] ;
          }
        }
        AA[row[i]][col [k]] = - AA[row[i]][col[k]] * AA[row[k]][col[k]] ;
      }
    } /* end inner reduction loop */
  } /* end main reduction loop */

  /* unscramble rows */
  for (j=0; j<n; j++)
  {
    for (i=0; i<n; i++)
    {
      temp[col[i]] = AA[row[i]][j];
    }
    for (i=0; i<n; i++)
    {
      AA[i][j] = temp[i] ;
    }
  }

  /* unscramble columns */
  for (i=0; i<n; i++)
  {
    for (j=0; j<n; j++)
    {
      temp[row[j]] = AA[i][col[j]] ;
    }
    for (j=0; j<n; j++)
    {
      AA[i][j] = temp[j] ;
    }
  }
  free(row);
  free(col);
  free(temp);
} /* end inverse */