#include <string.h>
#include "spnbox.h"
#include "../pnheaders/matrix.h"
#include "../pnheaders/pns.h"
#include "matrixmath.h"
#include "MemoryManager.h"

/*Subroutine declarations. See subroutine definitions for descriptions.*/
static int CheckParams(matrix* L, matrix* D, int** Tuc, int* TucCount, int** Tuo, int* TuoCount);
static void InitLPProblems(matrix *M1, matrix* M2, double** N1, double** N2, short** CType1, short** CType2, double** LoBounds, double** Cost, int* Tuc, int TucCount, int* Tuo, int TuoCount, matrix *D, MemoryManager * mem);
static ilpadm_r InitOther(matrix* CDeltas, matrix *L, matrix* D, int *Tuc, int TucCount, int* Tuo, int TuoCount, MemoryManager* mem);
static int TransformType(matrix* CDeltas, matrix* L, int* b, int* m0, int* Tuc, int TucCount, int* Tuo, int TuoCount, int Constraint);
static void BuildLP1(matrix *M1, matrix* CDeltas, double* Cost, int* Tuc, int TucCount, int* Tuo, int TuoCount, int Constraint);
static int RedoLP1(matrix* L, int* b, ipsolve_r* solution, int Constraint);
static void BuildLP1Alternate(matrix* M2, matrix* CDeltas, matrix* L, int *Tuc, int TucCount, int *Tuo, int TuoCount, int Constraint);
static void BuildLP2(matrix* M2, matrix* L, matrix* CDeltas, double* Cost, int* b, int* m0, int* Tuc, int TucCount, int* Tuo, int TuoCount, int Constraint);
static void ApplySolution(ilpadm_r* result, ipsolve_r* res, int Constraint, int Places, int Constraints);
static FinalizeSolution(ilpadm_r* result, matrix* L, int* b, int Places, int Constraints);
static void disp(char* beginning, char* parameter, char* end, int i, int p);

ilpadm_r ilpadm(matrix* L, int* b, matrix* D, int TucCount, int* Tuc, int TuoCount, int* Tuo, int* m0)
{
  /*The parameters to return.*/
  ilpadm_r result; 
  
  /*The return value from ipsolve*/
  ipsolve_r LPSolution;
  
  /*Variables. L is a constraint matrix.
  L has as many columns as there are places in the Petri net represented by D, 
	that is, as many columns as D has rows.
  L has as many rows as there are constraints.
  D is the incidence matrix of a Petri net. It has as many rows as there are 
	places and as many columns as there are transitions.
  CDeltas is a matrix of the change in constraint-met state (negative means more 
	met,	0 just met, positive means less met) that will occur when uncontrollable 
	or unobservable transitions fire. Each row holds all the deltas for a specific 
	constraint, and each column holds all the deltas for a specific transition*/ 
  int Constraints, Places, Transitions;  
  matrix CDeltas;
  
  /*Generic variables.*/
  int i = 0, j = 0;  
  
  /*The constraint matrices M1 & 2 , constraint vector N1 & 2, constraint type 
	vector CType1 & 2, and lower bounds LoBounds that will be used for the linear 
	programming problem when developing the transform. 2 is identical to 1 in all 
	cases except it has an extra constraint row. Which one we need is determined 
	per loop iteration.*/
  double *N1, *N2, *LoBounds, *Cost;
  short *CType1, *CType2;
  matrix M1, M2;
  
  /*Keep track of dynamically allocated memory so we can free it at the end.*/
  MemoryManager memory;  
  
  /*Initialize the return value to an error value.*/
  memset(&result, 0, sizeof(ilpadm_r));
  result.how = HOW_ERROR;  
  
  /*Check the parameters. Fail if any required parameters have not been given 
	and fill in defaults for any nonrequired parameters. L and D are required and 
	at least one of Tuc or Tuo is required. All other parameters are optional.*/
	if (! CheckParams(L, D, &Tuc, &TucCount, &Tuo, &TuoCount))
	{
		return result;
	}
	
  Constraints = NumberOfRows(*L);
  Places = NumberOfColumns(*L);
  Transitions = NumberOfColumns(*D);

  /*Create the memory manager.*/
  memory = CreateMemoryManager(10, 5, 10, 5);
  
  /*Allocate space for and initialize the descriptions of the two types of 
	linear programming problems that are possible when determining the
	transformation of any given constraint. Most of the problem will not change
	from constraint to constraint, so we can initialize within this function
	before the main constraint-processing loop.*/
	InitLPProblems(&M1, &M2, &N1, &N2, &CType1, &CType2, &LoBounds, &Cost, Tuc, TucCount,  Tuo, TuoCount, D, &memory);
	
	/*Allocate space for and initialize some other quantities that will be used
	mid-problem.*/
	result = InitOther(&CDeltas, L, D, Tuc, TucCount, Tuo, TuoCount, &memory);
  
  if(is_verbose() >= VRB_ILPADM)
  {
    printf("\n");
  }
  
  for (i = 0; i < Constraints; i++)
  {
    /*Status on each row defaults to OK.*/
    result.dhow[i] = HOW_OK;

		/*The ipsolve solution, used later, defaults to empty.*/
		memset(&LPSolution, 0, sizeof(ipsolve_r));
		
		/*See what kind of transformation needs to be done on the current constraint.*/		
		switch(TransformType(&CDeltas, L, b, m0, Tuc, TucCount, Tuo, TuoCount, i))
		{
		case -1:
			/*The constraint is infeasible. Set messages.*/
			result.dhow[i] = HOW_INFEASIBLE;
      disp("The transformation is infeasible due to the initial marking", "", ".", i, Constraints);			
			break;    
		case 0:
			/*No transformation is necessary.*/
			result.dhow[i] = HOW_OK;
			disp("No transformation is necessary", "", ".", i, Constraints);
			/*Copy the original marking constraint row, untransformed. Do so by 
			leaving this row of R1 all zeros and	making this element along R2's 
			diagonal equal to 1.*/
			SetMatrixEl(&result.R2, i, i, 1);
			break;		
		case 1:
			/*A type-1 transformation, one that ignores initial marking, should be done.
			Build the last few parameters of the linear programming problem.*/
			BuildLP1(&M1, &CDeltas, Cost, Tuc, TucCount, Tuo, TuoCount, i);
			/*Solve the linear programming problem.*/
			LPSolution = ipsolve(&M1, N1, Cost, 0, 0, LoBounds, CType1);			
			/*There are conditions under which we may need to resolve the problem with
			an additional constraint. Test for them.*/
			if (RedoLP1(L, b, &LPSolution, i))
			{
				DeallocateIpsolve(&LPSolution);
				/*Build the alternate linear programming problem and solve it.*/
				BuildLP1Alternate(&M2, &CDeltas, L, Tuc, TucCount, Tuo, TuoCount, i);
				LPSolution = ipsolve(&M2, N2, Cost, 0, 0, LoBounds, CType2);
				//Set the status message as appropriate.
				if (LPSolution.res)
				{
					result.dhow[i] = LPSolution.mhow;
				}
				else
				{
					result.dhow[i] = HOW_NOTSOLVED;
				}
			}
			/*Apply the resultes of the linear programming problem to the 
			transformation matrices.*/
			ApplySolution(&result, &LPSolution, i, Places, Constraints);			
			break;
		case 2:
			/*A type-2 transformation, one that takes into account the initial
			marking, should be done. Build the last few parameters of the linear
			programming problem.*/			
			BuildLP2(&M2, L, &CDeltas, Cost, b, m0, Tuc, TucCount, Tuo, TuoCount, i);
			/*Solve the LP.*/
			LPSolution = ipsolve(&M2, N2, Cost, 0, 0, LoBounds, CType2);
			result.dhow[i] = LPSolution.mhow;
			/*Apply the results to the transformation matrices.*/
			ApplySolution(&result, &LPSolution, i, Places, Constraints);
			break;
		}
  }
  
	/*Finalize the solution.*/
	FinalizeSolution(&result, L, b, Places, Constraints);  
	//Free memory and return.
  FreeMemory(&memory);
	return result;
}

/*******************************************************************************
CheckParams examines the parameters passed to ilpadm to make sure that required
parameters are present, does some basic validity checking such as making sure
matrix dimensions agree, and makes sure that all array pointers are null if 
array lengths are zero and vice-versa.*/
int CheckParams(matrix* L, matrix* D, int** Tuc, int* TucCount, int** Tuo, int* TuoCount)
{
  int Constraints, Places, Transitions;
	/*L and D are both required.*/
	if (! (L && D))
  {
    merror(0, "ilpadm: Required parameters were null pointers");
    return 0;
  }
  Constraints = NumberOfRows(*L);
  Places = NumberOfColumns(*L);
  Transitions = NumberOfColumns(*D);
  /*L must have the same number of columns as D does rows, since both of those 
	dimensions correspond to the number of places being constrained.*/
  if (Places != NumberOfRows(*D))
  {
    merror(0, "ilpadm: L does not have the same number of columns as D does rows");
    return 0;
  }
  /*Make sure that nullity of array pointers and zero in array lengths agree.*/
	if (! *Tuc)
	{
		*TucCount = 0;
	}
	else if (! *TucCount)
	{
		*Tuc = 0;
	}
	if (! *Tuo)
	{
		*TuoCount = 0;
	}
	else if (! *TuoCount)
	{
		*Tuo = 0;
	}
	/*At least one of Tuc or Tuo must be provided. If not, there is no 
	transformation to perform, since there are no unobservable or uncontrollable 
	transitions.*/
  if (!( *Tuc || *Tuo))
  {
    merror(0, "ilpadm: Both Tuc and Tuo are null pointers");
		return 0;
  }	
	return 1;
}

/*******************************************************************************
This function initializes the various parameters of the two alternate linear
programming problems, one of which will be solved for every constraint that has
to be transformed. The first group of parameters are pointers to the linear
programming parameters. The last 6 parameters are the uncontrollable transition
index list and index list length, the unraisable transition index list and index
list length, the incidence matrix, and the memory manager.
*/
void InitLPProblems(matrix *M1, matrix* M2, double** N1, double** N2, short** CType1, short** CType2, double** LoBounds, double** Cost, int* Tuc, int TucCount, int* Tuo, int TuoCount, matrix *D, MemoryManager * mem)
{
	int Places, i;
	Places = NumberOfRows(*D);
	
	/*Allocate space and initialize values.*/
  MAllocateMatrixType(mem, 2, M1, TucCount + TuoCount, Places + 1);
  MAllocateMatrixType(mem, 2, M2, TucCount + TuoCount + 1, Places + 1);
  *N1 = mcalloc(mem, (TucCount + TuoCount), sizeof(double));
  *N2 = mcalloc(mem, (TucCount + TuoCount + 1), sizeof(double));  
	*LoBounds = mcalloc(mem, (Places + 1), sizeof(double));
  *Cost = mcalloc(mem, (Places + 1), sizeof(double));
  *CType1 = mcalloc(mem, (TucCount + TuoCount), sizeof(short));
  *CType2 = mcalloc(mem, (TucCount + TuoCount + 1), sizeof(short));	

  /*If there are uncontrollable transitions, Tuc is an array of indexes to the 
	columns in D that contain them. Fill the first TucCount rows of M with these 
	columns. Since we are filling rows with columns, transpose M1 and M2 and 
	transpose them again after this block. Also, initialize the constraint types, 
	constraint vectors, and so forth.*/
	TransposeMatrix(M1);
	TransposeMatrix(M2);
  
	if(Tuc)
  {    
		for (i = 0; i < TucCount; i++)
    {
      /*Each integer problem constraint corresponding to an uncontrollable 
			transition is a <= type.*/
			(*CType1)[i] = -1;
      (*CType2)[i] = -1;

			/*Copy a single column, index Tuc[i], from D to the appropriate row of M1 
			and M2.*/
      CopyBlock(Places, 1, D, 0, Tuc[i], M1, 0, i);
			CopyBlock(Places, 1, D, 0, Tuc[i], M2, 0, i);			
		}
  }
  else
  {
    TucCount = 0;
  }
  /*Do the same thing for unobservable transitions.*/
  if(Tuo)
  {    
		for (i = 0; i < TuoCount; i++)
    {
			/*The constraint types corresponding to unobservable transitions are =. 
			They default to this. No need to change.*/
      /*Copy a single column, index Tuo[i], from D to the appropriate row of M1 
			and M2.*/
      CopyBlock(Places, 1, D, 0, Tuo[i], M1, 0, i + TucCount);
			CopyBlock(Places, 1, D, 0, Tuo[i], M2, 0, i + TucCount);			
    }
  }
  else
  {
    TuoCount = 0;
  }
  //The last constraint type 2 is another <=
  (*CType2)[TucCount + TuoCount] = -1;
  /*N1 is all zeros - no change necessary. N2 is all zeros except for the last 
	element.*/
  (*N2)[TucCount + TuoCount] = -1;
  /*LoBounds all zeros except the last one*/
  (*LoBounds)[Places] = 1;
  
	/*Un-transpose M1 and M2.*/
	TransposeMatrix(M1);
	TransposeMatrix(M2);
}

/*******************************************************************************
This function initializes parts of the return value ilpadm will use and various
intermediate variables used by the algorithm. CDeltas is a pointer to a matrix
with a row for each constraint and a column for each transition, each element of
which is the change in marking constraint satisfaction (negative being more
satisfied, positive being less) that occurs when the corresponding transition 
fires.*/
ilpadm_r InitOther(matrix* CDeltas, matrix *L, matrix* D, int *Tuc, int TucCount, int* Tuo, int TuoCount, MemoryManager* mem)
{ 
	int Constraints, Places, i, j;	
	ilpadm_r result;
	
	Constraints = NumberOfRows(*L);
	Places = NumberOfRows(*D);
	
	/*Initialize the return structure. Message defaults to OK, make space for 
	per-constraint messages.*/
	memset(&result, 0, sizeof(ilpadm_r));
	result.how = HOW_OK;
	result.dhow = tcalloc(Constraints, sizeof(char*));
	
	/*Allocate space for R1. It will be added to the constraint matrix, so it 
	should have a row for each constraint and a column for each place, just like 
	the original constraint matrix.*/
  AllocateMatrixType(2, &result.R1, Constraints, Places);
  /*Allocate space for R2. It will be a square matrix with the same number of 
	rows and columns as L has rows.*/
  AllocateMatrixType(2, &result.R2, Constraints, Constraints);
	/*We need a record of the changes to constraint-met state that occur for 
	each of the	uncontrollable and unobservable transitions. This corresponds 
	to the matrix	product of L * an incidence matrix containing only those 
	transitions. It has a row for each constraint and a column for each
	uncontrollable and unobervable transition. */
	MAllocateMatrixType(mem, 2, CDeltas, Constraints, TucCount + TuoCount);
	for (i = 0; i < Constraints; i++)
	{
		for (j = 0; j < TucCount; j++)
		{
			SetMatrixEl(CDeltas, i, j, MultiplyVector(L, D, i, Tuc[j], 0));
		}
		for (j = 0; j < TuoCount; j++)
		{
			SetMatrixEl(CDeltas, i, j + TucCount, MultiplyVector(L, D, i, Tuo[j], 0));
		}
	}
	return result;
}

/*******************************************************************************
TransformType examines the current constraint row and returns an integer
indicating what type of action should be taken. Return codes are:
-1	Transformation is infeasible.
 0	No transformation is necessary.
 1	A type-1 transformation, one ignoring the initial marking constraint, should
		be carried out.
 2	A type-2 transformation, one accounting for the initial marking constraint,
		should be carried out.
*/
int TransformType(matrix* CDeltas, matrix* L, int* b, int* m0, int* Tuc, int TucCount, int* Tuo, int TuoCount, int Constraint)
{
	int UCTDeltasNegative = 1, UOTDeltasZero = 1, i;

	/*First, check whether a transformation is feasible given the initial 
	marking and marking constraint vector (if they have been given)*/
	if (m0 && b && MultiplyVector(L, (matrix*) m0, Constraint, MULTIPLYVECTOR_ARRAY, 0) > b[Constraint])
	{
		return -1;
	}

/*UCTDeltasNegative if the current constraint is such that there are no 
	uncontrollable transitions the firing of which will ever cause it to change 
	from met to unmet. Set the UOTDeltasZero if the current constraint is such 
	that there are no unobservable transitions the firing of which will affect 
	the constraint-met state at all. Both of these flags get set if there are no
	appropriate transitions at all. Each column of the incidence matrix 
	corresponds to the change in marking that will occur if the corresponding 
	transition fires. So, a constraint row times an incidence column yields a 
	single number that will be negative if the transition firing causes the 
	constraint to become more met, positive if it causes it to become less met, 
	and 0 if it does not affect the state of the constraint. Each element of the 
	ConstraintDeltas matrix represents such a multiplication.	Rows are 
	constraints and columns are transitions.*/
	UCTDeltasNegative = 1;
	UOTDeltasZero = 1;    
	for (i = 0; i < TucCount; i++)
	{
		if(GetMatrixEl(CDeltas, Constraint, i) > 0)
		{
			UCTDeltasNegative = 0;
			break;
		}
	}
	for (i = 0; i < TuoCount; i++)	
	{
		if(GetMatrixEl(CDeltas, Constraint, i + TucCount))
		{
			UOTDeltasZero = 0;
			break;
		}
	}
	
	/*If the current constraint row is such that the firing of no 
	unobservable transition will change its met state and no uncontrollable 
	transition will ever cause it to become unmet, then no transformation is 
	necessary.*/
	if(UCTDeltasNegative && UOTDeltasZero)
	{
		return 0;
	}
	/*Otherwise, transformation is necessary. If there is an initial marking 
	vector, the constraint we are about to solve for will have to be met by 
	it. This constitutes an extra constrain row in the linear programming problem.
	Use the LP parameters with room for an extra row.*/
	else if (m0 && b)
	{
		return 2;
	}
	/*If not it we don't need the extra row. Use the smaller LP parameters.*/
	else
	{
		return 1;
	}
}

/*******************************************************************************
This function builds the final parts of the linear programming problem for a
type-1 transformation, those parts that vary from constraint row to constraint
row.
M1 and Cost are the constraint matrix and the cost vector the lp problem.
CDeltas, Tuc, TucCount, Tuo, and TuoCount are all the same as in ilpadm.
Constraint is the index of the current constraint row.
*/
void BuildLP1(matrix *M1, matrix* CDeltas, double* Cost, int* Tuc, int TucCount, int* Tuo, int TuoCount, int Constraint)
{
	int i, Places;
	Places = NumberOfColumns(*M1) - 1;
	/*Build the last column of the constraint vector.*/
	for (i = 0; i < TucCount; i++)
	{
		SetMatrixEl(M1, i, Places, GetMatrixEl(CDeltas, Constraint, i));
	}
	for (i = TucCount; i < TuoCount + TucCount; i++)
	{
		SetMatrixEl(M1, i, Places, GetMatrixEl(CDeltas, Constraint, i));
	}
	/*Fill the cost vector with ones.*/
	for (i = 0; i <= Places; i++)
	{
		Cost[i] = 1;
	}
}

/*******************************************************************************
This function examines the results of a type-1 transformation and returns a 1 if
the conditions under which the transformation needs to be redone with an 
additional constraint are met.
L is the constraint matrix, b is the constraint vector, solution is the
structure returned by the initial ipsolve call for the type-1 transformation,
and Constraint is the index of the current constraint row.*/
int RedoLP1(matrix* L, int* b, ipsolve_r* solution, int Constraint)
{
	/*We want to run a check. If b is nonempty and if the algorithm below 
	results in the setting of the flag, we will have to modify the problem 
	slightly and re-solve.*/
	int flag, i;
	if (solution->res && b)
	{
		flag = 1;
		for(i = 0; i < NumberOfColumns(*L); i++)
		{
			if (solution->res[i] + solution->res[NumberOfColumns(*L)] * GetMatrixEl(L, Constraint, i))
			{
				flag = 0;
				break;
			}
		}
		if (flag && b[Constraint] <= -1) return 1;	
	}
	return 0;
}

/*******************************************************************************
This function builds the alternate form of the linear programming problem
that is sometimes required for a type-1 transformation. Because this alternate
employs an extra constraint row in the lp problem, it actually uses those
lp parameters intended for a type-2 transformation since they all have an one 
more constraint row than those of a type-1 transformation.*/
void BuildLP1Alternate(matrix* M2, matrix* CDeltas, matrix* L, int *Tuc, int TucCount, int *Tuo, int TuoCount, int Constraint)
{
	int i, j, Places;
	/*We add one additional constraint to the linear programming problem. Its 
	coefficients are are a row of -1s followed by the inverse sum of all elements 
	in the current constraint vector row and its constant is a 1. Constraint type 
	should be set to greater than. N2 and ConstraintType2 are already set up this 
	way. We just need to fill the last row of M2 with the appropriate numbers and 
	put in the numbers we calculated for the last column of M2.*/
	j = 0;
	Places = NumberOfColumns(*M2) - 1;
	for(i = 0; i < Places; i++)
	{
		SetMatrixEl(M2, TucCount + TuoCount, i, 1);
		j += GetMatrixEl(L, Constraint, i);
	}
	SetMatrixEl(M2, TucCount + TuoCount, Places, j);
	/*Fill the last column of M2 as M1 was filled.*/
	for (i = 0; i < TucCount; i++)
	{
		SetMatrixEl(M2, i, Places, GetMatrixEl(CDeltas, Constraint, i));
	}
	for (i = 0; i < TuoCount; i++)
	{
		SetMatrixEl(M2, i + TucCount, Places, GetMatrixEl(CDeltas, Constraint, i + TucCount));
	}
}

/*******************************************************************************
This function builds the final parts of the linear programming problem for a
type-2 transformation, those parts that vary from constraint row to constraint
row.
M2 and Cost are the parameters of the linear programming problem.
Constraint is the index of the current constraint row.
The rest of the parameters have the same meaning as used in ilpadm.
*/
void BuildLP2(matrix* M2, matrix* L, matrix* CDeltas, double* Cost, int* b, int* m0, int* Tuc, int TucCount, int* Tuo, int TuoCount, int Constraint)
{
	int i, Places;
	Places = NumberOfColumns(*L);
	
	/*Build the last column of the constraint matrix (not counting the extra 
	row). Each element should be the change in current constraint 
	satisfaction that occurs when the transition corresponding to the 
	current row is fired. These values can be found in the products found 
	earlier.*/
	for (i = 0; i < TucCount; i++)
	{
		SetMatrixEl(M2, i, Places, GetMatrixEl(CDeltas, Constraint, i));
	}
	for (i = TucCount; i < TuoCount + TucCount; i++)
	{
		SetMatrixEl(M2, i, Places, GetMatrixEl(CDeltas, Constraint, i));
	}
	/*Build the cost vector and the last row of the constraint matrix. They 
	begin with the initial marking.*/
	for (i = 0; i < Places; i++)
	{
		/*The last row of the M2 constraint vector begins with the the initial 
		marking vector.*/
		SetMatrixEl(M2, TucCount + TuoCount, i, m0[i]);
		Cost[i] = m0[i];
	}
	/*The cost vector and the last constraint of the integer problem both 
	end with this number*/
	Cost[Places] = MultiplyVector(L, (matrix*) m0, Constraint, MULTIPLYVECTOR_ARRAY, 0) - b[Constraint] - 1;
	SetMatrixEl(M2, TucCount + TuoCount, Places, Cost[Places]);        
}

/*******************************************************************************
This function applies the results of an ipsolve call used to make a
transformation to the transformation matrices held in the ilpadm_r result 
structure.*/
void ApplySolution(ilpadm_r* result, ipsolve_r* LPSolution, int Constraint, int Places, int Constraints)
{
	int i;
	/*LPSolution.res contains results from type 1 or 2 transformation solutions 
	if it contains anything. Copy it into the appropriate places.*/
	if (LPSolution->res)
	{
		/*The last element of the solution gets copied to the current element in 
		the main diagonal of R2.*/
		SetMatrixEl(&result->R2, Constraint, Constraint, LPSolution->res[Places]);
		//The rest of the solution gets copied to the current row of R1.
		for(i = 0; i < Places; i++)
		{
			SetMatrixEl(&result->R1, Constraint, i, LPSolution->res[i]);
		}
		//Now, free the memory set aside by ipsolve for its result.
		free(LPSolution->res);
		LPSolution->res = 0;	
		
		/*If we are in verbose mode and the status is anything other than OK, we 
		want to show so.*/
		if(strcmp(LPSolution->mhow, HOW_OK))
		{
			disp("The problem is ", LPSolution->mhow, ".", Constraint, Constraints);
		}
	}
}

/*******************************************************************************
This function does the actual constraint transformation using the transformation
matrices built earlier. It also fills in the rest of the result structure
parameters.*/
FinalizeSolution(ilpadm_r* result, matrix* L, int* b, int Places, int Constraints)
{
	int i, j, x;
	/*Check for unsolved or infeasible solutions. The overall message should 
	prioritize to reporting infeasibilities, then report unsolved problems.*/
  for(i = 0; i < Places; i++)
  {
    if(result->dhow[i] == HOW_INFEASIBLE)
    {
      result->how = HOW_IMPOSSIBLE;
      break;
    }
    if(result->dhow[i] == HOW_NOTSOLVED)
    {
      result->how = HOW_NOTSOLVED;
    }
  }
  /*Get La according to formula: La = R1 + R2 * L*/
  /*Because R2 is a diagonal matrix, we can optimize things slightly here by 
	manually multiplying and examining only the diagonal, and then manually 
	adding.*/
	AllocateMatrixType(2, &result->La, Constraints, Places);
	for (i = 0; i < Constraints; i++)
	{
		x = 0;
		for(j = 0; j < Places; j++)
		{
			SetMatrixEl(&result->La, i, j, (GetMatrixEl(&result->R2, i, i) * GetMatrixEl(L, i, j)) + GetMatrixEl(&result->R1, i, j));
		}
	} 
  /*If b is available, get ba according to algorithm: ba = (R2 * a column vector 
	made of the elements of b + 1), each element then less 1. Do this with a loop 
	rather than matrix operations.*/
  if(b)
  {
    result->ba = tcalloc(Constraints, sizeof(int));
    for (i = 0; i < Constraints; i++)
    {
      result->ba[i] = (GetMatrixEl(&result->R2, i, i) * (b[i] + 1)) - 1;
    }
  }
}

/*******************************************************************************
disp displays debugging messages if the global is_verbose function indicates a
sufficiently high verbosity level. Otherwise it does nothing. Parameters are:
beginning, parameter, end - pointers to null-terminated strings that are
displayed on after another in the display string. Useful for adding beginning
and ending formats to ipsolve result message.
i - the current constraint row index.
p - the number of constraint rows. If this is 1 or less than the message will
not display constraint index information.*/
void disp(char* beginning, char* parameter, char* end, int i, int p)
{
  if(is_verbose() >= VRB_ILPADM)
  {
    if(p <= 1)
    {
      printf("ILPADM: %s%s%s\n", beginning, parameter, end);
    }
    else
    {
      printf("ILPADM, Constraint %d: %s%s%s\n", i, beginning, parameter, end);
    }
  }
}
