lpmex.c

matlab6 lpmex matlab6 lpmexma tlab6 lpmexmatlab6 lpmex

/* 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: [e-mail]
         [e-mail]

*/


#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 100

lprec	*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[])
{
	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]);
	}