/* compilation:  gcc pivot.c -o pivot */

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

double *ratios;
double **B_mat;
double EPSILON = 1e-5;

  FILE *d;
  char diary[40];

/***************************************************************/
void read_tableau(f_name, nb_rows, nb_cols, phase)
char *f_name;  
int *nb_rows; 
int *nb_cols; 
int *phase;
{
  int i, j;
  char buff[256];
  FILE *f;
  sprintf(diary,"%s.diary", f_name);

printf("f_name = %s\n",f_name);
  f = fopen(f_name, "r");
  d = fopen(diary, "w");

  fscanf(f,"%s", buff);

  while (strcmp(buff, "begin") != 0) {
    fscanf(f,"%s", buff);
  }
  fscanf(f, "%d %d", nb_rows, nb_cols);
 
  printf("#rows : %6d   #columns : %6d\n",*nb_rows, *nb_cols);
  fprintf(d,"#rows : %6d   #columns : %6d\n",*nb_rows, *nb_cols);
  fgets(buff, sizeof(buff),f);

  (*nb_rows)++;
  (*nb_cols)++;
  B_mat = (double**) malloc (sizeof(double*)*(*nb_rows));
  ratios = (double*) malloc (sizeof(double)*(*nb_rows));

  for (i=0; i<*nb_rows; i++)
    B_mat[i] = (double*) malloc (sizeof(double)*(*nb_cols));
        
  for (i=0; i<*nb_rows; i++) {
    
#ifdef TRACE
    printf("%d : ", i);
#endif

    for (j=1; j<*nb_cols; j++) {
      fscanf(f, "%lf", &(B_mat [i] [j])); 
        
#ifdef TRACE
      printf("%f  ", B_mat [i] [j]);
#endif
    }

#ifdef TRACE
    printf("\n");
#endif
  }
  fclose(f);

  *phase = 0;
} /* read_tableau */

/***************************************************************/
void
print_tableau(nb_rows, nb_cols, ph)
     int nb_rows;
     int nb_cols;
     int ph;
{ 
  int i, j;
  
  for(i=0; i<nb_rows; i++) 
  {
   for(j=0;j<nb_cols; j++)
   {
    if((B_mat[i][j] < EPSILON) && (B_mat[i][j] > -EPSILON)) 
    {
      B_mat[i][j] = 0;
    }
   }
  } 

  switch(ph) {
  
  case 0: 
    printf("PHASE 0\n");
    fprintf(d,"PHASE 0\n");
    break;

  case 1: 
    printf("PHASE 1\n");
    fprintf(d,"PHASE 1\n");
    break;

  case 2:
    printf("PHASE 2\n");
    fprintf(d,"PHASE 2\n");
    break;
  }
  printf("         ");
  fprintf(d,"         ");
  if(ph == 1) {
    printf("x 0      ");    
    fprintf(d,"x 0      ");    
  }
  for (j=1; j<nb_cols-1; j++) {
    printf("x%2d    ", j);
    fprintf(d,"x%2d    ", j);
  }
  printf("  RHS\n");
  fprintf(d,"  RHS\n");

  for (i=0; i<nb_rows-2;i++) {
    printf("%2d:   ", i+1);
    fprintf(d,"%2d:   ", i+1);

    if(ph == 1) {
      printf("%6.2f | ", B_mat[i][0]);      
      fprintf(d,"%6.2f | ", B_mat[i][0]);      
    }
    for (j=1; j<nb_cols-1; j++) {
      if((fabs(B_mat[i][j]) > EPSILON) && 
	 (fabs(B_mat[i][j]) < 0.005)) {
	printf("%5.1f* ", B_mat[i][j]);
	fprintf(d,"%5.1f* ", B_mat[i][j]);
      }
      else{
	printf("%6.2f ", B_mat[i][j]);
	fprintf(d,"%6.2f ", B_mat[i][j]);
      }
    }
    
    printf("| %6.2f\n", B_mat[i][nb_cols-1]);
    fprintf(d,"| %6.2f\n", B_mat[i][nb_cols-1]);
  }

  printf("      ");
  fprintf(d,"      ");
  if(ph == 1) {
    printf("--------");
    fprintf(d,"--------");
  }
  for (j=1; j<nb_cols; j++) {
    printf("--------");
    fprintf(d,"--------");
  }
  printf("\n");
  fprintf(d,"\n");
  printf("cost: ");
  fprintf(d,"cost: ");
  if(ph == 1) {
    printf("%6.2f | ", B_mat[nb_rows-2][0]);      
    fprintf(d,"%6.2f | ", B_mat[nb_rows-2][0]);      
  }
  for (j=1; j<nb_cols-1; j++) {
    if((fabs(B_mat[nb_rows-2][j]) > EPSILON) && 
       (fabs(B_mat[nb_rows-2][j]) < 0.005)) {
      printf("%5.1f* ", B_mat[nb_rows-2][j]);
      fprintf(d,"%5.1f* ", B_mat[nb_rows-2][j]);
    }
    else{
      printf("%6.2f ", B_mat[nb_rows-2][j]);
      fprintf(d,"%6.2f ", B_mat[nb_rows-2][j]);
    }
  }
  printf("| %6.2f", B_mat[nb_rows-2][nb_cols-1]);
  fprintf(d,"| %6.2f", B_mat[nb_rows-2][nb_cols-1]);
  printf("\n");
  fprintf(d,"\n");

  if(ph == 1) {
    printf("      ");
    fprintf(d,"      ");
    for (j=0; j<nb_cols; j++) {
      printf("--------");
      fprintf(d,"--------");
    }
    printf("\n");
    fprintf(d,"\n");
    printf("ph1 : ");
    fprintf(d,"ph1 : ");
    printf("%6.2f | ", B_mat[nb_rows-1][0]);
fprintf(d,"%6.2f | ", B_mat[nb_rows-1][0]);
    for (j=1; j<nb_cols-1; j++) {
      printf("%6.2f ", B_mat[nb_rows-1][j]);
fprintf(d,"%6.2f ", B_mat[nb_rows-1][j]);
    }
    printf("| %6.2f\n\n", B_mat[nb_rows-1][nb_cols-1]);
fprintf(d,"| %6.2f\n\n", B_mat[nb_rows-1][nb_cols-1]);
  }
  printf("\n");
fprintf(d,"\n");
} /* print_tableau */

/***************************************************************/
int main(int ac, char *av[]) {

  int i, j, nb_ineqs, nb_vars, r_piv = 1, c_piv = 1, error = 0, phase;
  char f_name_B[256];
  
  switch(ac)
  {
   case 1:
    printf("\nFile name for the tableau :\n");
    scanf("%s",f_name_B);
   break;
   case 2:
    sprintf(f_name_B,"%s",av[1]);
   break;
   default:
    fprintf(stderr,"Usage %s [File_name]\n",av[0]);
    exit(1);
   break;
  }

  read_tableau(f_name_B, &nb_ineqs, &nb_vars, &phase);

  while(((r_piv != 0) && (c_piv != 0)) || (phase < 3)) {

    print_tableau(nb_ineqs, nb_vars, phase);

    printf("Pivot (row# col#) ? (use 0 0 to change phase or end) : ");
fprintf(d,"Pivot (row# col#) ? (use 0 0 to change phase or end) : ");
    scanf("%d %d", &r_piv, &c_piv);
    fprintf(d,"%d %d", r_piv, c_piv);
    printf("\n");
fprintf(d,"\n");

    error = 0;
    if((r_piv < 0) || ((r_piv == 0) && (c_piv != 0)) || (r_piv > nb_ineqs-2)) {
      printf("ERROR: row index should be in the range [1, %d]\n", nb_ineqs-2);
fprintf(d,"ERROR: row index should be in the range [1, %d]\n", nb_ineqs-2);
      error = 1;
    } 

    if(phase == 1) {
      if((c_piv < 0) || (c_piv > nb_vars-2)) {
	printf("ERROR: column index should be in the range [0, %d]\n", 
	       nb_vars-2);
fprintf(d,"ERROR: column index should be in the range [0, %d]\n", 
	       nb_vars-2);
	error = 1;
      }
    }
    else {
      if((c_piv < 0) || (c_piv > nb_vars-2) || 
	 ((c_piv == 0) && (r_piv != 0))) {
	printf("ERROR: column index should be in the range [0, %d]\n", 
	       nb_vars-2);
fprintf(d,"ERROR: column index should be in the range [0, %d]\n", 
	       nb_vars-2);
	error = 1;
      }
    }
    if((c_piv == 0) && (r_piv == 0)) {
      phase++;
      error = 1;

      if(phase == 1) {
	for(i=0; i<nb_ineqs; i++) {
	  if(B_mat[i][nb_vars-1] < 0) {
	    B_mat[i][0] = -1;
	  }
	  else {
	    B_mat[i][0] = 0;
	  }
	}
	B_mat[nb_ineqs-1][0] = 1;
	for(j=1; j<nb_vars-1; j++) {
	  B_mat[nb_ineqs-1][j] = 0;
	}
	B_mat[nb_ineqs-1][nb_vars-1] = 0;
      }
    }
    else {
      if(!error && (B_mat[r_piv-1][c_piv] == 0)) {
	printf("ERROR: pivot can not be 0\n");
fprintf(d,"ERROR: pivot can not be 0\n");
	error = 1;
      }
    }
    if(!error) {
      for(j=0; j<nb_vars; j++) {
	if(c_piv != j) {
	  B_mat[r_piv-1][j] /= B_mat[r_piv-1][c_piv];
	}
      }
      B_mat[r_piv-1][c_piv] = 1;
      for(i=0; i<nb_ineqs; i++) {
	if(i != r_piv-1) {
	  for(j=0; j<nb_vars; j++) {
	    if(j != c_piv) {
	      B_mat[i][j] = B_mat[i][j]-B_mat[r_piv-1][j]*B_mat[i][c_piv];
	    }
	  }
	  B_mat[i][c_piv] = 0;
	}
      }
    }
  }
  return(0);
}