//======================================================= 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                                =
//=    7) Converts a Q matrix into a valid P matrix before solving          =
//=-------------------------------------------------------------------------=
//= 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: Ken 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)   =
//=           KJC (05/01/05) - Fixed a problem with solving Q matrices      =
//=           KJC (05/11/05) - Now can have multiple sequential comments    =
//=           KJC (06/16/06) - Cleaned-up some comments                     =
//===========================================================================
//----- 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
  int    col_num;                  // Column number for loop
  int    row;                      // Row number
  double val;                      // Value at row number
  double max_val;                  // Maximum sum value
  int    i;                        // Loop counter

  // 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
  max_val = 0;
  index = 0;
  do
  {
    scanf("%s", temp_string);

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

    // 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);

    // Find the maximum sum value (for Q matrix adjustment only)
    if (-Column[index].value > max_val)
      max_val = -Column[index].value;

    // Increment index
    index++;
  } while (1);

  // If Q type then need to convert to a valid P matrix
  if (Matrix_type == Q_TYPE)
  {
    index = 0;
    for(col_num=0; col_num<Size; col_num++)
    {
      // Loop for each column
      do
      {
        row = Column[index].row_num;
        val = Column[index].value;

        // Compute the adjusted value
        val = val / max_val;
        if (col_num == row)
          val = val + 1.0;

        // Save the adjusted value
        Column[index].value = val;

        // Increment index
        index = index + 1;

        // if row is -1 then no more values in this column
        if (row < 0) break;
      } 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 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
      sum = sum + (Pi[row] * val);

    } while (1);

    // Determine next value of pi_temp
    pi_temp[col_num] = sum;
  }

  // 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);
}