//======================================================= file = iter.c =====
//=  Solver for P*pi = pi and Q*pi = 0 for sparse P and Q                   =
//=   - Uses iterative methods                                              =
//===========================================================================
//=  Notes:                                                                 =
//=    1) Input from input file "in.dat" to stdin (see example below)       =
//=        * Note special format of input file, "-1" denotes end of         =
//=          column, "-999" denotes last entry, and "&" bounds comments     =
//=    2) Reads matrix Type ('P' or 'Q') as first input value and then      =
//=       reads matrix Size (0...(Size - 1)) as second input value.  No     =
//=       initial comment is allowed.                                       =
//=    3) Output is to stdout                                               =
//=    4) The matrix is NOT checked for validity (i.e., that for a P        =
//=       matrix every row sums to 1.0 and for a Q matrix that every row    =
//=       sums to 0.0).  A non-valid matrix will cause unpredicatable       =
//=       results.                                                          =
//=    5) Assumes that sum of pi is always 1.0.                             =
//=    6) Uses direct iteration for P matrix and Gauss-Seidel for Q matrix. =
//=-------------------------------------------------------------------------=
//= Example "in.dat" file:                                                  =
//=                                                                         =
//=   P 3                                                                   =
//=   & Example from Kleinrock, Volume 1, pages 31 - 33                     =
//=     ---------------- column #0    &                                     =
//=   1       0.25                                                          =
//=   2       0.25                                                          =
//=   -1      zero                                                          =
//=   & ---------------- column #1    &                                     =
//=   0       0.75                                                          =
//=   2       0.25                                                          =
//=   -1      zero                                                          =
//=   & ---------------- column #2    &                                     =
//=   0       0.25                                                          =
//=   1       0.75                                                          =
//=   2       0.50                                                          =
//=   -1      zero                                                          =
//=   & ---------------- END FLAG     &                                     =
//=   -999    zero                                                          =
//=-------------------------------------------------------------------------=
//= Example output (for above "in.dat"):                                    =
//=                                                                         =
//=   -------------------------------------------------- iter.c -----       =
//=     Iteration count = 10 -- Convergence mean = 0.333333                 =
//=     Iteration count = 20 -- Convergence mean = 1.66151e-06              =
//=     Iteration count = 30 -- Convergence mean = 3.20145e-12              =
//=     Iteration count = 40 -- Convergence mean = 0                        =
//=   ---------------------------------------------------------------       =
//=     Results for P matrix - size = 3 and tolerance = 1e-15               =
//=       Pi[0] = 0.2                                                       =
//=       Pi[1] = 0.28                                                      =
//=       Pi[2] = 0.52                                                      =
//=   ---------------------------------------------------------------       =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 iter.c, cl iter.c, gcc iter.c -lm                         =
//=-------------------------------------------------------------------------=
//=  Execute: iter < in.dat                                                 =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (02/23/02) - Genesis (from "original" iter.c from 1991)   =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <stdio.h>                 // Needed for printf()
#include <math.h>                  // Needed for fabs()
#include <stdlib.h>                // Needed for exit(), atoi(), and atof()
#include <string.h>                // Needed for strcmp()

//----- Defines -------------------------------------------------------------
#define P_TYPE           0         // P matrix type is 0
#define Q_TYPE           1         // Q matrix type is 1

#define MAX_SIZE        50         // Max size of Q array (0...(MAX_SIZE - 1))
#define MAX_ENTRY     1000         // Max number of non-zero entries
#define TOLERANCE  1.0e-15         // Converence tolerance for pi values
#define NUM_LOOP        10         // Iterations between convergence checks

//----- Globals -------------------------------------------------------------
struct Col_type                    // Data structure for matrix columns
{
  int    row_num;                  // ** Row number for this column entry
  double value;                    // ** Value for this column entry
};
struct Col_type Column[MAX_ENTRY]; // Columns of matrix
int             Matrix_type;       // Matrix type (P or Q)
int             Size;              // Size of matrix  (0...(Size - 1))
int             Count;             // Count for iterations
int             Done;              // Flag for convergence
double          Pi[MAX_SIZE];      // pi vector
double          Pi_old[MAX_SIZE];  // pi vector from previous iteration
double          Converge_sum;      // Convergence sum

//----- Prototypes ----------------------------------------------------------
void load_matrix(void);            // Load matrix
void iterate(void);                // Perform one iteration
int  check_convergence(void);      // Check for convergence of Pi to Pi_old

//===========================================================================
//=  Main program                                                           =
//===========================================================================
void main(void)
{
  int     i;                       // Loop counter

  // Load the matrix
  printf("-------------------------------------------------- iter.c -----\n");
  load_matrix();

  // Initialize Pi
  for (i=0; i<Size; i++)
    Pi[i] = 1.0 / Size;

  // Main loop for iterative solution
  Count = Done = 0;
  do
  {
    Count = Count + 1;

    for (i=0; i<NUM_LOOP; i++)
      iterate();

    Done = check_convergence();

    printf("  Iteration count = %d -- Convergence mean = %g \n",
      (Count * NUM_LOOP), (Converge_sum / Size));

  } while (Done == 0);

  // Output results
  printf("---------------------------------------------------------------\n");
  if (Matrix_type == P_TYPE)
    printf("  Results for P matrix - ");
  else
    printf("  Results for Q matrix - ");
  printf("size = %d and tolerance = %g \n", Size, TOLERANCE);
  for (i=0; i<Size; i++)
    printf("    Pi[%d] = %g \n", i, Pi[i]);
  printf("---------------------------------------------------------------\n");
}

//===========================================================================
//=  Function to load matrix from stdin                                     =
//===========================================================================
void load_matrix(void)
{
  char    temp_string[256];        // Temporary string variable
  int     index;                   // Index value

  // Read first entry, must be matrix type
  scanf("%s", temp_string);
  if ((temp_string[0] == 'P') || (temp_string[0] == 'p'))
    Matrix_type = P_TYPE;
  else if ((temp_string[0] == 'Q') || (temp_string[0] == 'q'))
    Matrix_type = Q_TYPE;
  else
  {
    printf("*** ERROR - Matrix type ('P' or 'Q') is not first entry \n");
    exit(1);
  }

  // Read second entry, must be matrix size
  scanf("%s", temp_string);
  Size = atoi(temp_string);
  if ((Size <=0) || (Size > MAX_SIZE))
  {
    printf("*** ERROR - Size is invalid (Size = %d) \n", Size);
    exit(1);
  }

  // Main loop to read contents of file into matrix
  index = 0;
  do
  {
    scanf("%s", temp_string);

    // This handles a comment
    if (strcmp(temp_string, "&") == 0)
    {
      do
      {
        scanf("%s", temp_string);
      } while (strcmp(temp_string, "&") != 0);
      scanf("%s", temp_string);
    }

    // Input row number
    Column[index].row_num = atoi(temp_string);
    if (Column[index].row_num == -999) break;

    // Input value for this row number
    scanf("%s", temp_string);
    Column[index].value = atof(temp_string);

    // Increment index
    index++;

  } while (1);
}

//===========================================================================
//=  Function to perform one iteration                                      =
//===========================================================================
void iterate(void)
{
  int    col_num;                  // Column number for loop
  int    row;                      // Row number
  double val;                      // Value at row number
  double diag_val;                 // Current diagonal value
  double sum;                      // Sum variable
  double pi_temp[MAX_SIZE];        // Temporary pi vector
  int    i;                        // Loop counter

  // Loop to perform one interation of Gauss-Seidel
  i = 0;
  for(col_num=0; col_num<Size; col_num++)
  {

    // Loop for each column
    sum = 0.0;
    do
    {
      row  = Column[i].row_num;
      val  = Column[i].value;

      i = i + 1;
      if (i >= MAX_ENTRY)
      {
        printf("*** ERROR - exceeded MAX_ENTRY elements (MAX_ENTRY = %d) \n",
          MAX_ENTRY);
        exit(1);
      }

      // row = -1 means no more values in this column
      if (row < 0) break;

      // An iteration
      if (Matrix_type == P_TYPE)
        sum = sum + (Pi[row] * val);
      else
      {
        if (row != col_num)
          sum = sum + (Pi[row] * val);
        else
          diag_val = val;
      }

    } while (1);

    // Determine next value of pi_temp
    if (Matrix_type == P_TYPE)
      pi_temp[col_num] = sum;
    else
      pi_temp[col_num] = (-sum) / diag_val;
  }

  // Normalize pi_temp
  sum = 0;
  for (i=0; i<Size; i++)
    sum = sum + pi_temp[i];
  for (i=0; i<Size; i++)
    pi_temp[i] = pi_temp[i] / sum;

  // Copy pi_temp into Pi
  for(i=0; i<Size; i++)
    Pi[i] = pi_temp[i];
}

//===========================================================================
//=  Function to check convergence of Pi to Pi_old                          =
//===========================================================================
int check_convergence(void)
{
  int    i;                        // Loop counter

  // Loop for convergence sum
  Converge_sum = 0.0;
  for (i=0; i<Size; i++)
    Converge_sum = Converge_sum + fabs(Pi[i] - Pi_old[i]);

  // Copy current Pi into Pi_old
  for (i=0; i<Size; i++)
    Pi_old[i] = Pi[i];

  // Return a "1" if converged
  if (Converge_sum <= TOLERANCE)
    return(1);
  else
    return(0);
}