/************************************************************************
* retroLP solves linear programs using the simplex method.             *
* Copyright (C) 2001, 2003, 2005 Gavriel Yarmish & Richard Van Slyke.  *
*                                                                      *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or    *
* (at your option) any later version.                                  *
*                                                                      *
* This program is distributed in the hope that it will be useful,      *
* but WITHOUT ANY WARRANTY; without even the implied warranty of       *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
* GNU General Public License for more details.                         *
*                                                                      *
* You should have received a copy of the GNU General Public License    *  
* along with this program; if not, write to the Free Software          *
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,           *
*                                                                      *
************************************************************************/


// This is the driver for the simplex method.
// CHANGES
// 4/10/05  Put in facility for user modification of parameters.
// 3/233/05 Integrated readmps in here
// 3/9/05   Defined a class LongClock for extending the clock
//          to longer runs, also put in clock time. Finally,
//          if reinverting at end of Phase I restarted counter
//          for next reinversion. Took out code for measuring
//          clock time. It was insignificant.
// 9/9/02   Integrate pivot and column choice
// 8/4/01   Added save_state to dump to .lp_ file for output
//          It is called "simplex infile" where infile is in the
//          DPLP binary format and describes the program instance.
// 7/19/02  Reversed matrix storage to column major
// 7/11/00  Added reinversion capability
// 6/27/00  Replaced runtime test for column choice
//          rule with compiler directive CCR
// 6/26/00  Replaced realtime debugging switches
//          with compiler directive DEBUG

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <errno.h>
#include <time.h>
#include "simplex.h"
#include "long_clock.h"
#include <string>
#include <string.h>
#include <algorithm>
#include <cctype>
#include <iostream>
#include <fstream>
#include <sstream>
#include "params.h"


LP_state * readmps(const char *, FILE *);
int to_file(FILE *, LP_state *);
int to_state(FILE *, short, LP_state **);
void free_a(double **, long, long);
short check_ph1(LP_state * );
void save_state(int, char *, LP_state *);
inline int add_to_ph1(LP_state *,long );
inline int sub_from_ph1(LP_state *,long );
void tableau_out(LP_state *);
long choose_column(LP_state *, double *,short * , double);
long choose_row(LP_state *, long , short,double *, iter_result *, double);
long change_state(short , LP_state *, long ,long , double, iter_result,
  double *, short *, double);
void report(long, LP_state *);
long refresh(FILE *, LP_state *, short *);
char *out_type(var_type);
char *out_range(var_range);
char *out_result(iter_result);
void out_istate(long,LP_state *);
void out_jstate(long,LP_state *);
void free_arrays(LP_state *);
FILE * cget_file(char *, int);
FILE * cbget_file(const char *, int);
void empty_file(const char *, FILE * );
int  get_file(std::ifstream&, const std::string&, int);
int UpDate_p(std::string , std::string );
void p_setup();
int p_parse(std::ifstream & );
void ip_parse();
void ShortHelp();
void LongHelp();
int KeyWordHelp(int, int);


// START MAIN HERE
// Main returns 0 if optimum achieved; -1 if unbounded; +1 if
// infeasible; -2 if maximum number of iterations exhausted

int main(int argc, char *argv[])
{

  // Analyze command line arguments

  // 5 Cases:
  //    1. argv[i] starts with a letter; then it is
  //       an input file name.
  //    2. argv[i] is "--"; we switch to interactive
  //       input.
  //    3. argv[i] starts with - but is not -- or
  //       -pfile. We get the value from argv.
  //    4. argv[i] is -pfile filename. We get values
  //       from the file filename
  // Examples:
  // "retroLP file.mps" runs the problem defined by the
  //   mps file using default parameters.
  // "retroLP --" takes parameters interactively from 
  //   cin (stdin).
  // "retroLP -pfile pfile.txt" takes parameter changes
  //   from the file pfile.txt
  // "retroLP -iter_max 10000 file.mps" changes the 
  //   maximum number of iterations to 10000 and then
  //   runs file.mps.
  ///////////////////////////////////////////////////
  // WE ASSUME ALL FILE NAMES START WITH A LETTER!

  std::string s, sold;
  bool expect_parameter_file = FALSE;
  bool expect_parameter_arg = FALSE;
  char c;
  int e;
  std::ifstream ph;

  if(argc <= 1)
  {
    std::cout << "\nretroLP requires an argument; use retroLP -help for help" << std::endl;
    return 11;
  }

  // Load parameter defaults
  p_setup();

  for(int i=1; i < argc; i++)
  {
    sold = s;
    s = argv[i];
    if(expect_parameter_arg)
      {
	// This should be the value part of a parameter change
	expect_parameter_arg = FALSE;
	// Negative return -> recycle s (e.g., for HELP)
	if(UpDate_p(sold,s) >= 0)continue;
      }
    // Pick out first character
    c = s[0];
    // Test for file name (mps or parameter)
        if(isalpha(c))
        {
	  // If expected parameter starts with an alpha
	  // Could be mps file or parameter file
	  if(expect_parameter_file == TRUE)
	    {
	      expect_parameter_file = FALSE;	      
	      // Open parameter file
	      std::cout << "Parameter File: " << s <<std::endl;
	      // Open file with handle ph
	      e = get_file(ph, s, 1);
	      // Loop through lines until eof.
	      e = p_parse(ph);
	      if(e)
		{
		  std::cerr << "Couldn't read parameter file: " << s << std::endl;
		  exit(e);
		}
	      // Each time call function to change parameter
	      // (same function to be used throughout)
	    }
	  else
	    {
	      //Otherwise assume it is a mps file to be run with current parameters
	      std::cout << "Input file: " << s << std::endl;
	      goto start_run;
	    }
        }
	
	// Test for intitial dash
	if(c == '-')
	  {
	    // Gobble up the dash
	    // If starts with dash, we can convert to upper case
		  transform(s.begin(),s.end(),s.begin(), toupper);
	    s.erase(s.begin());
	    // Then process
	    c = s[0];
	    // Look for "--"
	    if(c == '-')
	      // Go to interactive mode
	      {
		std::cout << "retroLP version 0.5, Copyright (C) 2001. 2003, "<<std::endl
			  << "2005 Gavriel Yarmish and Richard Van Slyke. "<<std::endl
			  << "retroLP comes with ABSOLUTELY NO WARRANTY; for "<<std::endl
			  << "details type `show w'. This is free software, "<<std::endl
			  << "and you are welcome to redistribute it under "<<std::endl
			  << "certain conditions; type `show c' for details."<<std::endl
			  << std::endl;
		std::cout << "Interactive mode--Enter parameter name (no dash), value and enter:" 
			  << std::endl;
		std::cout << "Control d or \"quit\" to stop interactive mode." << std::endl;
		ip_parse();
		continue;
	      }
	    // Look for file mode
	    if(strcmp(s.c_str(), "PFILE")==0)      
	      {
		expect_parameter_file = TRUE;
		std::cout << "File mode: " << s << std::endl;
		continue;
	      }
	    // Parameter change from argument
	    // Pick up value next iteration
	    expect_parameter_arg = TRUE;
	  }
  }
  // If we run out of arguments or end with  an action with no arguments
  if(expect_parameter_file||expect_parameter_arg)
    {
      if(strcmp(s.c_str(), "HELP")==0)
	{
	  ShortHelp();
	  exit(0);
	}
      std::cerr << "Open parameter change or missing parameter file name!" << std::endl;
      exit(2);
    }
  // Missing .mps file
  exit(3);
 start_run:
  LP_state *lp;
  short phase, old_phase=0, up_flag, r1=0, r2=0, error_spot=0;
  iter_result nonbasicflag;
  int r;

  long  ip, jp, jn = 0, frefresh;
  double dx, valj, delta = del0;
  clock_t clock0;
  clock0=clock();
  time_t time0;
  time0 = time(NULL);
  printf("\n%s \n", ctime(&time0));
  
  // OPEN FILE (name now comes in string "s")
  FILE *  bfile;
  // Open up a binary file to hold LP structure on disk
  bfile = cbget_file("bfile.lpb",13);
  if(bfile == NULL)return(13);

  // to build the dynamic memory structure, 
  // and copy it to a binary file.
  // The dynamic memory structure is freed at the end of main.
  // It is pointed to by lp.
  
  
  // See if minimizing or maximizing
  if(max)maxfactor = -1.0;else maxfactor = 1.0;

  // Open mps file
  lp = readmps(s.c_str(), bfile);

  // If you construct your own LP_state, you can start here

  if(lp == NULL)
    {
      std::cerr << "\nError in input file: " << s << std::endl;
      error_spot = 14;
      goto clean;
    }
  
  // -------------- START SIMPLEX METHOD ------------
  // Check for Phase 1
  phase=check_ph1(lp);

  // Maximum number of iterations before reinverting basis
  // The setting of frefresh is a result of tradeoffs
  // For numerical accuracy you want it to be fairly
  // small; to avoid cycling you want it large.
  // RESET_FREQ represents the cycling interest
  // and lp->m*refresh_factor the numerical accuracy interest.
  // This is our compromise
  frefresh=( reset_freq < lp->m*refresh_factor )?reset_freq
    :lp->m*refresh_factor;
  frefresh=(1000>frefresh)?1000:frefresh;
  printf("Reinversion frequency is %ld\n",frefresh);
  switch(ccr)
    {
    case 1: printf("We are using the \"classical\" column choice rule\n");
      break;
    case 2: printf("We are using the \"steepest edge\" column choice rule\n");
      break;
    case 3: printf("We are using the \"greatest change\" column choice rule\n");
      break;
    }
  if(l_start > 0)
    {
      printf("======================================\n");
      printf(" AT START: \n");
    }       
  report(l_start,lp);
  
  // Reinitialize long clocks
  tlong.resetLongClock();
  ptlong.resetLongClock();
  pclock = 0;
  
  
  // MAIN ITERATION STARTS HERE
  
  for(iter=1; iter<= iter_max; iter++)
    {
      // Check to see if LongClock needs updating
      tlong.reinitLongClock();
      // Update count for reinversion
      fiter++;
      // If enough iterations reinvert basis and reset delta
      if(refresh_off == 0)
	{
	  if(fiter%frefresh == 0)
	    {
	      old_phase=phase;
	      printf("======================================\n");
	      printf("RESET called at iteration %ld, major %ld, re %ld\n",
		     iter, miter, riter);
	      printf("Time = %s\n", tlong.getInterval());
	      refresh(bfile, lp, &phase);
	      fiter = 1;
	      delta = del0;
	      jn = 0;
	      if(phase!=old_phase)
		printf("Warning: phase changed from %ld to %ld during reset\n",
		       old_phase, phase);
	    }
	}
      delta+=tau;
      switch(phase)
	{
	case 1:
	  // --------------------- PHASE I -------------------
	  
	  // Pass delta for greatest change column choice
	  // If jn is non-zero then we got value from pivot+ and don't
	  // need to look again.
	  if(jn == 0)jp = choose_column(lp, &valj, &up_flag, delta );
	  else jp = jn;
	  // Check for end of Phase I
	  if(jp==0) // No column offers improvement
	    {
	      // Reset once
	      if(refresh_off == 0)
		{
		  if(!r1)
		    {
		      r1++;
		      printf("======================================\n");
		      printf("RESET at end of Ph 1: iteration %ld, major %ld, re %ld\n",
			     iter, miter, riter);
		      printf("Time = %s\n", tlong.getInterval());
		      refresh(bfile, lp, &phase);
		      fiter = 1;
		      delta = del0;
		      jn = 0;
		      jp = choose_column(lp, &valj, &up_flag, delta );
		      if((jp!=0)&&(phase==1))
			{
			  printf(" Reverting to Phase I\n");
			  break;
			}
		    }
		}
	      else
		{
		  if(!r1)
		    {
		      r1++;
		      jn = 0;
		      phase = check_ph1(lp);
		      jp = choose_column(lp, &valj, &up_flag, delta );
		      if((jp!=0)&&(phase==1))
			{
			  printf(" Reverting to Phase I without refresh\n");
			  break;
			}
		    }
		}
	      if(lp->b[lp->m+1]>-epsb) // Make sure feasible
		{
		  phase=2;
		  lp->mm = 0;
		  if(l_ph1>0)
		    {
		      printf("\n======================================\n");
		      printf("END OF PHASE 1 at iteration %ld, major %ld, re %ld\n",
			     iter, miter, riter);
		      report(l_ph1,lp);
		      jn = 0;
		    }
		}
	      else // Infeasible
		{
		  printf("\n======================================\n");
		  printf("INFEASIBLE SOLUTION at iteration %ld, major %ld, re %ld\n",
			 iter, miter, riter);
		  report(l_end, lp);
		  r = 1;
		  goto clean;
		}
	      
	      break;
	    }
	  ip = choose_row(lp, jp, up_flag, &dx, &nonbasicflag, delta);
	  // We shouldn't get ip==-1; if so, refresh. If still -1 abort.
	  if(ip==-1)
	    {
	      printf("Error in phase 1; unbounded column %ld at iteration %ld\n",jp, iter);
	      error_spot = 222;
	      goto clean;
	    }
	  if((iter%k_iter == 0) && (l_iter >= 2))
	    {
	      printf("Phase %hd,  %s, row: %ld, column: %ld, ",
		     phase, up_flag?"increase":"decrease", ip, jp);
	      printf("dx: %10.2e, valj: %10.2e, \nresult: %s\n",
		     dx, valj, out_result(nonbasicflag));
	      printf("pivot element: %12.8e\n", lp->a[jp][ip]);
	      printf("BEFORE\n");
	      out_istate(ip, lp);
	      out_jstate(jp, lp);
	    }
	  assert((ip <= lp->m) && (ip >= 0));
	  assert((jp <= lp->n) && (jp > 0));
	  jn = change_state(phase, lp, ip, jp, dx, nonbasicflag, &valj,
			    &up_flag, delta);
	  if((iter%k_iter == 0) && (l_iter >= 2))
	    {
	      printf("AFTER\n");
	      out_istate(ip, lp);
	      out_jstate(jp, lp);
	    }
	  break;
	case 2:
	  // --------------------- PHASE II ---------------------
	  // Pass delta for greatest change column choice
	  // If jn is non-zero then we got value from pivot+ and don't
	  // need to look again.
	  if(jn == 0) jp = choose_column(lp, &valj, &up_flag, delta );
	  else jp = jn;
	  // Check for optimum value
	  
	  if(jp==0) // If jp is still 0 no column offers improvement
	    {
	      if(!refresh_off)
		// Refresh once
		{
		  if(!r2)
		    {
		      printf("======================================\n");
		      printf("RESET called at end of Ph2, iteration %ld, major %ld, re %ld\n",
			     iter, miter, riter);
		      printf("Time = %s\n", tlong.getInterval());
		      refresh(bfile, lp, &phase);
		      fiter = 1;
		      r2++;
		      if(phase==1)
			{
			  printf(" False end of phase 2 at iteration %ld, major %ld, re %ld\n",
				 iter, miter, riter);
			  jn = 0;
			  break;
			}
		      else
			{
			  // try column choice again
			  jp = choose_column(lp, &valj, &up_flag, delta);
			  if(jp!=0)
			    {
			      printf(" False end of phase 2 at iteration %ld, major %ld, re %ld\n",
				     iter, miter, riter);
			      printf(" reverted to phase 2 in refresh\n");
			      jn = 0;
			      break;
			    }
			  
			}
		    }
		}
	      if(l_ph1>0)
		{
		  report(l_ph1,lp);
		}
	      printf("\n======================================\n");
	      printf("OPTIMAL SOLUTION at iteration %ld, major %ld, re %ld\n",
		     iter, miter, riter);
	      report(l_end, lp);
	      r = 0;
	      goto clean;
	    }
	  ip = choose_row(lp, jp, up_flag, &dx, &nonbasicflag, delta);
	  // Check for unbounded solution
	  if(ip < 0)
	    {
	      printf("\n======================================\n");
	      printf("UNBOUNDED SOLUTION at iteration %ld, major %ld, re %ld, Column %ld \n",
		     iter, miter, riter, -jp);
	      report(l_end,lp);
	      r = -1;
	      goto clean;
	    }
	  if((iter%k_iter == 0) && (l_iter >= 2))
	    {
	      printf("Phase %hd,  %s, row: %ld, column: %ld, ",
		     phase, up_flag?"increase":"decrease", ip, jp);
	      printf("dx: %10.2e, valj: %10.2e, \nresult: %s\n",
		     dx, valj, out_result(nonbasicflag));
	      printf("pivot element: %12.8e\n", lp->a[jp][ip]);
	      printf("BEFORE\n");
	      out_istate(ip, lp);
	      out_jstate(jp, lp);
	    }
	  assert((ip <= lp->m) && (ip >= 0));
	  assert((jp <= lp->n) && (jp > 0));
	  jn = change_state(phase, lp, ip, jp, dx, nonbasicflag, &valj,
			    &up_flag, delta);
	  if((iter%k_iter == 0) && (l_iter >= 2))
	    {
	      printf("AFTER\n");
	      out_istate(ip, lp);
	      out_jstate(jp, lp);
	    }
	  break;
	default:
	  printf("\n wrong value %d for phase\n",phase);
	  error_spot = 20;
	  goto clean;
	} // end switch
      if(iter%k_iter == 0)report(l_iter,lp);
    } // end for
  printf("\n======================================\n");
  printf("Reached MAXIMUM NUMBER OF ITERATIONS\n");
  r=-2;
  report(l_end, lp);
  // ------------- CLEAN UP, RELEASE FILES AND MEMORY ----------
  
 clean:
  
  if(lp != NULL)
    {
      printf("Wall clock time: %15.3f\n", difftime(time(NULL),time0));
      printf("Running time: %15.3f\n",
	     ((double)(clock()-clock0))/CLOCKS_PER_SEC);
      // Output of pivot timing
      printf("Running time using long clock: %s\n", tlong.getInterval());
      printf("Pivot time: %15.3f\n",((double)pclock/CLOCKS_PER_SEC));
      printf("Pivot time using long clock: %s\n", ptlong.getInterval(pclock));
      printf("+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=\n");
      // Free arrays
      free_arrays(lp);
    }
  empty_file("bfile.lpb", bfile);
  printf("\n<><><><><><><><><><><><><><><><><><><>\n");
  if(error_spot == 0)exit(r);
  else exit(-error_spot);
}
//===============================CHECK_PH1=======================
short check_ph1(LP_state *s)
{
  long i, j;
  short p=2;
  // Zero out Phase I objective
  for (j=0; j<=s->n; j++)s->a[j][s->m + 1]=0.0;
  s->b[s->m + 1] = 0.0;
  
  // LOOK AT NON-BASICS
  
  // Fill in nrange based on x and bounds
  for (j=1; j<=s->n; j++)
    {
      switch(s->ntype[j])
	{
	case FREE:
	  s->nrange[j]= BETWEEN;
	  break;
	  
	case LBOUNDED:
	  if(s->x[j] < s->nl[j]-epsb)
	    {
	      s->nrange[j]= BELOW;
	      p=1;
	      s->a[j][s->m+1]+=1.0;
	    }
	  else
	    {
	      if((s->x[j] >= s->nl[j]-epsb)&&(s->x[j]<=s->nl[j]+epsb))
		s->nrange[j]= ATLOWER;
	      else s->nrange[j]= BETWEEN;
	    }
	  break;
	  
	case UBOUNDED:
	  if(s->x[j] > s->nu[j]+epsb)
	    {
	      s->nrange[j]= ABOVE;
	      p=1;
	      s->a[j][s->m+1]+=-1.0;
	    }
	  else
	    {
	      if((s->x[j] >= s->nu[j]-epsb)&&(s->x[j]<=s->nu[j]+epsb))
		s->nrange[j]= ATUPPER;
	      else s->nrange[j]= BETWEEN;
	    }
	  break;
	  
	case BOTHBOUNDS:
	  if(s->x[j] < s->nl[j]-epsb)
	    {
	      s->nrange[j]= BELOW;
	      p=1;
	      s->a[j][s->m+1]+=1.0;
	    }
	  else
	    {
	      if((s->x[j] >= s->nl[j]-epsb)&&(s->x[j]<=s->nl[j]+epsb))
		s->nrange[j]= ATLOWER;
	      else if (s->x[j] < s->nu[j]-epsb)s->nrange[j]= BETWEEN;
	      else if (s->x[j] > s->nu[j]+epsb)
		{
		  s->nrange[j]= ABOVE;
		  p=1;
		  s->a[j][s->m+1]+=-1.0;
		}
	      else s->nrange[j] = ATUPPER;
	    }
	  break;
	  
	case FIXED:
	  if(s->x[j] < s->nl[j]-epsb)
	    {
	      s->nrange[j]= BELOW;
	      p=1;
	      s->a[j][s->m+1]+=1.0;
	    }
	  else
	    {
	      if((s->x[j] >= s->nl[j]-epsb)&&(s->x[j]<=s->nl[j]+epsb))
		s->nrange[j]= ATBOTH;
	      else
		{
		  s->nrange[j]= ABOVE;
		  p=1;
		  s->a[j][s->m+1]+=-1.0;
		}
	    }
	  break;
	  
	default:
	  printf("\n Error in brange: j = %ld\n",j);
	  exit(501);
	}
    }
  // LOOK AT BASICS
  // Fill in brange based on b and bounds
  // We don't allow BASIC variables to be AT bounds
  // except ATBOTH (11/28/99)
  
  
  for (i=1; i<=s->m; i++)
    {
      switch(s->btype[i])
	{
	case FREE:
	  s->brange[i]= BETWEEN;
	  break;
	  
	case LBOUNDED:
	  if(s->b[i] < s->bl[i]-epsb)
	    {
	      s->brange[i]= BELOW;
	      p=1;
	      add_to_ph1(s, i);
	    }
	  else
	    {
	      s->brange[i]= BETWEEN;
	    }
	  break;
	  
	case UBOUNDED:
	  
	  if(s->b[i] > s->bu[i]+epsb)
	    {
	      s->brange[i]= ABOVE;
	      p=1;
	      sub_from_ph1(s, i);
	    }
	  else
	    {
	      s->brange[i]= BETWEEN;
	    }
	  break;
	  
	case BOTHBOUNDS:
	  if(s->b[i] < s->bl[i]-epsb)
	    {
	      s->brange[i]= BELOW;
	      p=1;
	      add_to_ph1(s, i);
	    }
	  else
	    {
	      if(s->b[i] > s->bu[i]+epsb)
		{
		  s->brange[i]= ABOVE;
		  p=1;
		  sub_from_ph1(s, i);
		}
	      else
		{
		  s->brange[i]=BETWEEN;
		}
	    }
	  break;
	  
	case FIXED:
	  if(s->b[i] < s->bl[i]-epsb)
	    {
	      s->brange[i]= BELOW;
	      p=1;
	      add_to_ph1(s, i);
	    }
	  else
	    {
	      if((s->b[i] >= s->bl[i]-epsb)&&(s->b[i]<=s->bl[i]+epsb))
		s->brange[i]= ATBOTH;
	      else
		{
		  s->brange[i]= ABOVE;
		  p=1;
		  sub_from_ph1(s, i);
		}
	    }
	  break;
	  
	default:
	  printf("\n Error in brange: i = %ld\n",i);
	  exit(502);
	}
    }
  
  if(p==2)s->mm=0;
  else s->mm=s->m+1;
  return p;
}
//=========================ADD_TO_PH1===========================
inline int add_to_ph1(LP_state *s,long i)
{
  long j;
  for(j=0; j<=s->n; j++) s->a[j][s->m + 1]+=s->a[j][i];
  s->b[s->m + 1]+=s->b[i];
  return 0;
  
}
//=========================SUB_FROM_PH1===========================
inline int sub_from_ph1(LP_state *s,long i)
{
  long j;
  for(j=0; j<=s->n; j++) s->a[j][s->m + 1]-=s->a[j][i];
  s->b[s->m + 1]-=s->b[i];
  return 0;
  
}
//========================SAVE_STATE==============================
// This routine is used to save an LP_state structure to a file
// It is an extension of an .lpb file, with additional information
// for output, checkpointing, etc.
void save_state(int r, char *filename, LP_state * lp)
{
  FILE* file;

  file = fopen(filename, "wb");
  if(file == NULL)
    {
      printf("\nError opening file for save: %s\n", filename);
      return;
    }
  
  to_file(file, lp);
  
  // Write out additional information here
  // Termination state
  // r = 0 if optimum achieved; -1 if unbounded; +1 if
  // infeasible; -2 if maximum number of iterations exhausted
  fwrite(&r, sizeof(int), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  
  // Variables for output (definition)
  fwrite(&l_start, sizeof(short), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&l_ph1, sizeof(short), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&l_end, sizeof(short), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&l_iter, sizeof(short), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&k_iter, sizeof(short), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  //Variables for statistics
  fwrite(&iter, sizeof(long), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&iter, sizeof(long), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&miter, sizeof(long), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&riter, sizeof(long), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&cntnfx, sizeof(long), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&cntcc, sizeof(long), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  // Support for measuring pivot times
  fwrite(&pclock, sizeof(clock_t), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&pstart, sizeof(clock_t), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  fwrite(&pstop, sizeof(clock_t), 1, file);
  if(feof(file)+ferror(file)>0) goto error;
  return;
  
 error:
  printf("Error in save_state -- feof, ferror: %d %d \n",
         feof(file), ferror(file));
  clearerr(file);
  exit(199);
}

//========================= P_PARSE ===========================

int p_parse(std::ifstream & pin)
{
  int r = 1;
  std::string line, code, s38d(38,'-');
  std::string  value;
  while(getline(pin, line))
    {
      code = "";
      value ="";
      // Read lines
      std::cout << s38d << std::endl;
      std::cout << line << std::endl;
      // Parse lines
      std::istringstream stream(line);
      stream >> code;
	  transform(code.begin(),code.end(),code.begin(), toupper);
      stream >> value;
      // Call change parameter function
      if(!(code == ""))UpDate_p(code, value);
    }
  if(pin.eof())r=0;
  return r;
}

//========================= IP_PARSE ===========================

// Handles interactive input of parameters

void ip_parse()
{
  std::string code, value, line;
  std::istringstream stream(line);
  // User enters a line here
  while(getline(std::cin, line))
    {
      // Parse lines
      std::istringstream stream(line);
      code = "";
      stream >> code;
	  transform(code.begin(),code.end(),code.begin(), toupper);
      // Check for quit
      if(strcmp(code.c_str(), "QUIT")==0) break;
      value = "";
      stream >> value;
      std::cout << "Line is: " << line << std::endl;
      std::cout << "Code: " << code << "; value: " << value << std::endl;
      // Return on eof or quit
      UpDate_p(code, value);
    }
  return;
}

