/* lpmex.c  Mex file interface to the lp_solve 2.3 toolkit. Please see  README_MEX for more information.  This file and the matlab interface code is Copyright (c) 1995,  Jeffrey C. Kantor. All Rights Reserved.    Email: jeff@control.cheg.nd.edu         cchen@darwin.helios.nd.edu*/#include <string.h>#include "mex.h"/* mex.h defines REAL as 0 to use as a flag, and the lp toolkit   defines REAL as a double. Therefore we need to undef REAL *//*#undef REAL*/#include "lpkit.h"/* Declare a global lp record */#define LPMAX 100lprec	*lp[LPMAX];int	lp_last,lp_nr;int	result;/* Some working storage */short	initialized=FALSE;char	*cmd; /* command */char	*errmsg; /* error message *//* An exit function to clean up the lp structure.   Registered with matlab and called by matlab on exit */void ExitFcn(){	int	i;	for (i=0; i<=lp_last; i++) {	  if (lp[i] != NULL)	          delete_lp(lp[i]);	}	free(cmd);	free(errmsg);}/* Function to get a real scalar with error checking */REAL GetRealScalar(	const mxArray	*pm){	if ((mxGetM(pm) == 1) && (mxGetN(pm)==1)	  && (mxIsNumeric(pm)) && (!mxIsComplex(pm)) ) {		return(mxGetScalar(pm));	} else {		mexErrMsgTxt("Expecting a scalar argument.");	}}/* Function to get len elements from a matlab vector. Matrix   can be either full or sparse. Elements are stored in indices   1..n  Errors out if the matlab vector is not length len */void GetRealVector(	const mxArray	*pm,	REAL	*vec,	int	len){	int	j,k,m,n;	int	*ir, *jc;	double	*pr;	m = mxGetM(pm);	n = mxGetN(pm);	if (   ((m == 1) && (len != n)) || ((n == 1) && (len != m))	  || !mxIsNumeric(pm) || mxIsComplex(pm)) {		mexErrMsgTxt("invalid vector function.");	}	pr = mxGetPr(pm);	if (!mxIsSparse(pm)) {		for (k=0; k<len; k++) {			vec[k+1] = pr[k];		}	} else if (mxIsSparse(pm)) {		jc = mxGetJc(pm);		ir = mxGetIr(pm);		for (j=0; j<n; j++) {			for (k=jc[j]; k<jc[j+1]; k++) {				vec[1+ir[k]+j*m] = pr[k];			}		}	} else {		mexErrMsgTxt("Can't figure out this matrix.");	}}/* lp_handle validation test */int handle_valid(        int handle){        if (lp[handle] == NULL || handle < 0 || handle > lp_last)	   return(0);        else           return(1);}/* mexFunction is the Matlab Interface */void mexFunction(        int	nlhs,        mxArray	*plhs[],        int	nrhs,        const mxArray	*prhs[]){        /* Start Change*/        mxArray *pa;        /* End Change*/	int	h,i,j,k,l,n,m;	int	row,col;	short	type;	REAL	value;	REAL	*yout;	REAL	*vec;	int	*ir, *jc;	double	*pr;	FILE    *fp;	nstring filename;	char    *str, *str1;	if (!initialized) {		/* Register the Exit Function */		if (0 != mexAtExit(ExitFcn)) {			mexErrMsgTxt("Failed to register exit function.\n");		}        	initialized = TRUE;		/* Allocate a string array to store command */		cmd = calloc(NAMELEN,sizeof(char));		/* Allocate a string array to store error message */		errmsg = calloc(200,sizeof(char));		/* Initialize the lp array, and pointer to the last lp */		for (i = 0; i < LPMAX; i++) 			lp[i] = NULL;		lp_last = -1;	}	/* Get the first argument as a string matrix */	if (nrhs < 1) {		mexErrMsgTxt("At least one command is required.");	}	mxGetString(prhs[0], cmd, mxGetN(prhs[0])+1);	/* Now call the required part of the lp toolkit *//* functions without using lp_handle */	/* [lp_handle] = lpmex('make_lp', rows, cols) */	if (strcmp("make_lp",cmd)==0) {		if (nrhs != 3) {			mexErrMsgTxt("make_lp requires 2 arguments.");		}		row = GetRealScalar(prhs[1]);		col = GetRealScalar(prhs[2]);		for (i=0; i<=lp_last; i++) {		  if (lp[i] == NULL)		      break;		}		if (i <= lp_last) {		  lp[i] = make_lp(row,col);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = i;		}		else if (lp_last+1 >= LPMAX ) {			mexErrMsgTxt("Cannot allocate any more lps.");		}		else {		  lp_last++;		  lp[lp_last] = make_lp(row,col);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = lp_last;		}	}	/* [lp_handle] = lpmex('read_lp_file', filename, verbose, lp_name) */	/* Make lp from a LP fromat file */	else if (strcmp("read_lp_file",cmd)==0) {		if (nrhs < 2) {			mexErrMsgTxt("read_lp_file requires at least 1 argument.");		}		mxGetString(prhs[1], filename, NAMELEN);		if ((fp = fopen(filename, "r")) == NULL) {		   mexErrMsgTxt("read_lp_file can't read file.");		}		type = (nrhs == 3) ? 1 : 0;		if (nrhs < 4) {		   strcpy(filename, "unnamed");		} 		else {		   mxGetString(prhs[3], filename, NAMELEN);		}		for (i=0; i<=lp_last; i++) {		  if (lp[i] == NULL)		      break;		}		if (i <= lp_last) {		  lp[i] = read_lp_file(fp, type, filename);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = i;		}		else if (lp_last+1 >= LPMAX ) {			mexErrMsgTxt("Cannot allocate any more lps.");		}		else {		  lp_last++;		  lp[lp_last] = read_lp_file(fp, type, filename);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = lp_last;		}		fclose(fp);	}		/* [lp_handle] = lpmex('read_mps', filename, verbose) */	/* Make lp from a MPS fromat file */	else if (strcmp("read_mps",cmd)==0) {		if (nrhs < 2) {			mexErrMsgTxt("read_mps requires 1 or 2 arguments.");		}		mxGetString(prhs[1], filename, NAMELEN);		if ((fp = fopen(filename, "r")) == NULL) {		   mexErrMsgTxt("read_mps can't read file.");		}		type = (nrhs == 3) ? 1 : 0;		for (i=0; i<=lp_last; i++) {		  if (lp[i] == NULL)		      break;		}		if (i <= lp_last) {		  lp[i] = read_mps(fp, type);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = i;		}		else if (lp_last+1 >= LPMAX ) {			mexErrMsgTxt("Cannot allocate any more lps.");		}		else {		  lp_last++;		  lp[lp_last] = read_mps(fp, type);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = lp_last;		}		fclose(fp);	}	/* [handle_vec] = lpmex('print_handle') */	else if (strcmp("print_handle",cmd)==0) {		j=0;	        for (i=0; i <= lp_last; i++) {		  if (lp[i] != NULL) {		     j++;		  }                }		if (j > 0) {		  plhs[0] = mxCreateDoubleMatrix(j,1,0);		  pr = mxGetPr(plhs[0]);		  k=0;		  for (i=0; i<= lp_last; i++) {		    if (lp[i] != NULL) {                     pr[k] = i;		     k++;		    }		  }		}		else {		 mexErrMsgTxt("Cannot find existing lp_handle."); 		}	} 	/* lpmex('demo') *//* This demo is not a necessary part of the program	else if (strcmp("demo",cmd)==0) {		lp_last++;		lp[lp_last] = make_lp(0,4);		str_add_constraint(lp[lp_last],"3 2 2 1",LE,4);		str_add_constraint(lp[lp_last],"0 4 3 1",GE,3);		str_set_obj_fn(lp[lp_last],"2 3 -2 3");		print_lp(lp[lp_last]);		solve(lp[lp_last]);		print_solution(lp[lp_last]);		plhs[0] = mxCreateDoubleMatrix(1,1,0);		pr = mxGetPr(plhs[0]);		pr[0] = lp_last;	} *//* functions with at least one argument */	else if (nrhs < 2 ) {                strcpy(errmsg,cmd);                strncat(errmsg,": Unimplemented or requires at least 1 argument.",180);                mexErrMsgTxt(errmsg);	} 	/* getting lp_handle */	else if (  !( (h=GetRealScalar(prhs[1])) + 1 )   )     {        strcpy(errmsg,cmd);        strncat(errmsg,": lp handle can not be -1.",180);	        		mexErrMsgTxt(errmsg);	} 	/* lp_handle validation */	else if ( !handle_valid(h) )     {        strcpy(errmsg,cmd);        strncat(errmsg,": Invalid lp handle.",180);	        		mexErrMsgTxt(errmsg);	} 	/* lpmex('auto_scale', lp_handle) */	else if (strcmp("auto_scale",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("auto_scale requires 1 argument.");		auto_scale(lp[h]);	}	/* [lp_handle] = lpmex('copy_lp', lp_handle) */		else if (strcmp("copy_lp", cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("copy_lp requires 1 argument.");		if (lp_last+1 >= LPMAX ) {			mexErrMsgTxt("Cannot allocate any more lps.");		}		for (i=0; i<=lp_last; i++) {		  if (lp[i] == NULL)		      break;		}		if (i <= lp_last) {		  lp[i] = copy_lp(lp[h]);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = i;		}		else if (lp_last+1 >= LPMAX ) {			mexErrMsgTxt("Cannot allocate any more lps.");		}		else {		  lp_last++;		  lp[lp_last] = copy_lp(lp[h]);		  plhs[0] = mxCreateDoubleMatrix(1,1,0);		  pr = mxGetPr(plhs[0]);		  pr[0] = lp_last;		}	}	/* lpmex('delete_lp', lp_handle) */	else if (strcmp("delete_lp", cmd)==0) {		if (nrhs != 2) {			mexErrMsgTxt("delete_lp requires 1 argument.");		}		delete_lp(lp[h]);                lp[h] = NULL;	}	/* [basis] = lpmex('get_basis', lp_handle) */	else if (strcmp("get_basis",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("get_basis requires 1 argument.");		if (lp[h]->basis_valid) {			plhs[0] = mxCreateDoubleMatrix(lp[h]->columns,1,0);			pr = mxGetPr(plhs[0]);			for (i=0; i<lp[h]->columns; i++) {				pr[i] = lp[h]->basis[i+1];			}		}	}	/* [cost_vec]=lpmex('get_reduced_costs',lp_handle) */ 	else if (strcmp("get_reduced_costs",cmd)==0) {	      if (nrhs != 2) {		 mexErrMsgTxt("get_reduced_costs needs 1 argument1.");	      }	      plhs[0] = mxCreateDoubleMatrix(lp[h]->sum,1,0);	      pr = mxGetPr(plhs[0]);	      get_reduced_costs(lp[h], pr);	}	/* [obj,x,duals] = lpmex('get_solution', lp_handle) */	else if (strcmp("get_solution",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("get_solution requires 1 argument.");		if (result == 0) {			plhs[0] = mxCreateDoubleMatrix(1,1,0);			pr = mxGetPr(plhs[0]);			pr[0] = lp[h]->best_solution[0];			if (nlhs > 1) {				plhs[1] = mxCreateDoubleMatrix(lp[h]->columns,1,0);				pr = mxGetPr(plhs[1]);				for (i=0; i<lp[h]->columns; i++) {					pr[i] = lp[h]->best_solution[lp[h]->rows+i+1];				}			}			if (nlhs > 2) {				plhs[2] = mxCreateDoubleMatrix(lp[h]->rows,1,0);				pr = mxGetPr(plhs[2]);				for (i=0; i<lp[h]->rows; i++) {					pr[i] = lp[h]->duals[i+1];				}			}		}	}	/* lpmex('print_duals', lp_handle) */	else if (strcmp("print_duals",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("print_duals requires 1 argument.");		print_duals(lp[h]);	}	/* lpmex('print_lp', lp_handle) */	else if (strcmp("print_lp",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("print_lp requires 1 argument.");		print_lp(lp[h]);	}	/* lpmex('print_scales', lp_handle) */	else if (strcmp("print_scales",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("print_scales requires 1 argument.");		print_scales(lp[h]);	}	/* lpmex('print_solution', lp_handle) */	else if (strcmp("print_solution",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("print_solution requires 1 argument.");		print_solution(lp[h]);	}	/* lpmex('reset_basis', lp_handle) */	else if (strcmp("reset_basis",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("reset_basis requires 1 argument.");		reset_basis(lp[h]);	}	/* lpmex('set_maxim', lp_handle) */	else if (strcmp("set_maxim",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("set_maxim requires 1 argument.");		set_maxim(lp[h]);	}	/* lpmex('set_minim', lp_handle) */	else if (strcmp("set_minim",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("set_minim requires 1 argument.");		set_minim(lp[h]);	}	/* [result] = lpmex('solve', lp_handle) */	else if (strcmp("solve",cmd)==0) {		if (nrhs != 2)			mexErrMsgTxt("solve requires 1 argument.");		result = solve(lp[h]);		plhs[0] = mxCreateDoubleMatrix(1,1,0);		pr = mxGetPr(plhs[0]);		pr[0] = result;		if (result == OPTIMAL) {		  if (lp[h]->verbose == TRUE) {		    print_solution(lp[h]);		    printf("Branch & Bound depth: %d\n",lp[h]->max_level);		    printf("Nodes processed: %d\n",lp[h]->total_nodes);		    printf("Simplex pivots: %d\n",lp[h]->total_iter);          }		}		if (result == INFEASIBLE)		  printf("This problem is infeasible\n");		if (result == UNBOUNDED)		  printf("This problem is unbounded\n");		if (result == FAILURE)		  printf("lp_solve failed\n");	}	/* lpmex('unscale', lp_handle) */