glpkmex.c

This routine calls the glpk library to solve a LP/MIP problem. A typical LP problem has following s

/*----------------------------------------------------------------------
 *
 * Copyright (C) 2001-2005, Nicolo' Giorgetti, All rights reserved.
 * E-mail: <[e-mail]>.
 *
 * This file is part of GLPK (GNU Linear Programming Kit).
 *
 * GLPK 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, or (at your option)
 * any later version.
 *
 * This part of code is distributed with the FURTHER condition that it 
 * can be compiled and linked with the Matlab libraries and it can be 
 * used within the Matlab environment.
 *
 * GLPK 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 GLPK; see the file COPYING. If not, write to the Free
 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 * 02111-1307, USA.
 *
 *-----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <string.h>
#include <setjmp.h>
#include "mex.h"
#include "glpk.h"

#include "glpksets.h"
#include "glpkfun.h"

#define	SENSE_IN   prhs[0]
#define	C_IN	    prhs[1]
#define	A_IN	    prhs[2]
#define	B_IN	    prhs[3]
#define  CTYPE_IN   prhs[4]
#define  LB_IN	    prhs[5]
#define  UB_IN      prhs[6]
#define  VARTYPE_IN prhs[7]
#define  PARAM      prhs[8]
#define  SOLVER_IN  prhs[9]
#define  SAVE_IN    prhs[10]


/* Output Arguments */
#define	 XMIN_OUT     plhs[0]
#define	 FMIN_OUT     plhs[1]
#define	 STATUS_OUT   plhs[2]
#define   EXTRA_OUT    plhs[3]

/* This constant is used to discriminate not null elements 
   of the constant matrix A. For instance, if A[i,j] > TOL
   A[i,j] is considered a not-null element of the matrix, 
   otherwise it is assumed NULL.
*/ 
#define TOL 1e-14


void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
   int    sense;
   double *tmp=NULL;
   double *c=NULL;
   double *A=NULL;
   double *b=NULL;
   double *bcopy=NULL; 
   double *lb=NULL;
   double *ub=NULL;
   char *vartype=NULL;
   int mrowsc=0;
   int mrowsA=0;
   char *ctype=NULL;
   int save_pb=0;
   int lpsolver=1;
   int error;
   double *xmin=NULL;
   double *fmin=NULL;
   double *status=NULL;
   double *lambda=NULL;
   double *redcosts=NULL;
   double *time=NULL;
   double *mem=NULL;
   char errmsg[1024];
   int *freeLB=0;
   int *freeUB=0;
   int nfields=0;
   const char **extranames;
   mxArray *mxlambda, *mxredcosts, *mxtime,*mxmem;
   double *rdtmp=NULL;
   mxArray *mxtmp;
   int *vartype2=NULL;
   int jmpret;
   int i;
   int size;
   int flagLB=0; /* flag used to free LB if allocated with mxCalloc */
   int flagUB=0; /* flag used to free UB if allocated with mxCalloc */
   

   /* row and column sets of non-zero constraint coefficients */
   int nz;  /* number of non-zero coefficients */
   int *rn=NULL;
   int *cn=NULL;
   double *a=NULL;  /* non-zero coefficients */

   /*
      flag to identify the type of problem:
             isMIP=0 <=> LP problem
             isMIP=1 <=> MIP problem
   */
   int isMIP = 0;

   if(nrhs < 1){
      mexPrintf(header);
      mexPrintf(version);
      mexPrintf(copyright);
      mexPrintf(syntax);
      return;
   }
   if(nrhs < 4) mexErrMsgTxt("At least 4 inputs required (SENSE,C,A,b)");
   if(nlhs < 2) mexErrMsgTxt("2 outputs required");


   /* 1st Input. Sense of optimization. */
   if (!mxIsNumeric(SENSE_IN)
       || (mxGetNumberOfDimensions(SENSE_IN) > 2)
       || (mxGetM(SENSE_IN) != 1)
       || (mxGetN(SENSE_IN) != 1)
       || mxIsComplex(SENSE_IN)
       || ((tmp = mxGetPr(SENSE_IN)) == NULL)
       || ((*tmp != 1) && (*tmp != -1))
      )
       mexErrMsgTxt("SENSE must be either 1 or -1.");
   else
       sense = (int) (*tmp);

   /* 2nd Input. A column array containing the objective function
                 coefficients.
   */
   if (!mxIsNumeric(C_IN)
	    || (mxGetNumberOfDimensions(C_IN) > 2)
	    || ((mrowsc = mxGetM(C_IN)) < 1)
	    || (mxGetN(C_IN) != 1)
	    || mxIsComplex(C_IN)
	    || ((c= mxGetPr(C_IN)) == NULL)
		)
		mexErrMsgTxt("C must be a real valued column vector.");

   /* 3rd Input. A matrix containing the constraints coefficients. */
	if (!mxIsNumeric(A_IN)
	    || (mxGetNumberOfDimensions(A_IN) > 2)
	    || ((mrowsA = mxGetM(A_IN)) < 1)
	    || (mxGetN(A_IN) != mrowsc)
	    || mxIsComplex(A_IN)
	    || ((A = mxGetPr(A_IN)) == NULL)
		) {
		sprintf(errmsg,"A must be a real valued %d by %d matrix.",mrowsA, mrowsc);
		mexErrMsgTxt(errmsg);
	}else{

	  if(!mxIsSparse(A_IN)){
	    int i,j;
	   /* int nrcount; */

	    rn=(int *)mxCalloc(mrowsA*mrowsc+1,sizeof(int));
	    cn=(int *)mxCalloc(mrowsA*mrowsc+1,sizeof(int));
	    a=(double *)mxCalloc(mrowsA*mrowsc+1,sizeof(double));

		/*
		nz=0;
	    for(i=0;i<mrowsA;i++){
	      for(j=0;j<mrowsc;j++){
		      if(A[i+j*mrowsA]!=0) nz++;
	      }
	    }
		rn=(int *)mxCalloc(mrowsA*mrowsc+1,sizeof(int));
	    cn=(int *)mxCalloc(nz+1,sizeof(int));
	    a=(double *)mxCalloc(nz+1,sizeof(double)); */

	    
	    nz=0;
	    for(i=0;i<mrowsA;i++){
	      for(j=0;j<mrowsc;j++){
		      if(fabs(A[i+j*mrowsA])> TOL){
                    nz++;
		            rn[nz]=i+1;
		            cn[nz]=j+1;
		            a[nz]=A[i+j*mrowsA];
		      }
	      }
	    }
	  }else{
	    int i,j;
	    int *jc,*ir;
	    double *pr;
	    int nelc,count,row;

	    /* NOTE: nnz is the actual number of nonzeros and is stored as the
       last element of the jc array where the size of the jc array is the
       number of columns + 1 */
	    nz = *(mxGetJc(A_IN) + mrowsc);
	    jc = mxGetJc(A_IN);
	    ir = mxGetIr(A_IN);
	    pr = mxGetPr(A_IN);

       rn=(int *)mxCalloc(nz+1,sizeof(int));
	    cn=(int *)mxCalloc(nz+1,sizeof(int));
	    a=(double *)mxCalloc(nz+1,sizeof(double));

       count=0; row=0;
	    for(i=1;i<=mrowsc;i++){
	      nelc=jc[i]-jc[i-1];
	      for(j=0;j<nelc;j++){
		      count++;
		      rn[count]=ir[row]+1;
		      cn[count]=i;
		      a[count]=pr[row];
		      row++;
	      }
	    }
	  }

	}

   /* 4th Input. A column array containing the right-hand side value
	         for each constraint in the constraint matrix.
   */
	if (!mxIsNumeric(B_IN)
	    || (mxGetNumberOfDimensions(B_IN) > 2)
	    || (mxGetM(B_IN) != mrowsA)
	    || (mxGetN(B_IN) != 1)
	    || mxIsComplex(B_IN)
	    || ((b = mxGetPr(B_IN)) == NULL)
		) {
		sprintf(errmsg,"B must be a real valued %d by 1 column vector.",mrowsA);
		mexErrMsgTxt(errmsg);
	}

   /* 5th Input. A column array containing the sense of each constraint
                 in the constraint matrix.
   */
  if ((nrhs > 4) && (mxGetM(CTYPE_IN) != 0) && (mxGetN(CTYPE_IN) != 0)) {
  if (!mxIsChar(CTYPE_IN)
       || (mxGetNumberOfDimensions(CTYPE_IN) > 2)
       || (mxGetM(CTYPE_IN) != mrowsA)
       || (mxGetN(CTYPE_IN) != 1)
       || mxIsComplex(CTYPE_IN)
      ){
	sprintf(errmsg,"CTYPE must be a char valued %d by 1 column vector.",mrowsA);
	mexErrMsgTxt(errmsg);
       } else {
	 /* int i,size; */

	size = mxGetNumberOfElements(CTYPE_IN) + 1;

	/* Allocate enough memory to hold the converted string. */
	ctype = mxCalloc(size, sizeof (char));

	/* Copy the string data from string_array_ptr and place it into buf. */
	if (mxGetString(CTYPE_IN, ctype, size) != 0)
	  mexErrMsgTxt("Could not convert string data.");

	/* checking if the input is made only of  F, U, S, L and D */
	for (i = 0; i < size - 1; i++) {
	  if (     (ctype[i] != 'F') && (ctype[i] != 'U')
		&& (ctype[i] != 'S') && (ctype[i] != 'L')
		&& (ctype[i] != 'D'))
	          mexErrMsgTxt("CTYPE must contain only F,U,S,L and D");
	}
       }
  }

   /* 6th Input. An array of at least length numcols containing the lower
            	  bound on each of the variables.
   */
    if ((nrhs > 5) && (mxGetM(LB_IN) != 0) && (mxGetN(LB_IN) != 0)) {
    if (!mxIsNumeric(LB_IN)
	|| (mxGetNumberOfDimensions(LB_IN) > 2)
	|| (mrowsc != mxGetM(LB_IN))
	|| (mxGetN(LB_IN) != 1)
	|| mxIsComplex(LB_IN)
	|| ((lb = mxGetPr(LB_IN)) == NULL)
       ) {
      sprintf(errmsg,"LB must be a real valued %d by 1 column vector.",mrowsc);
      mexErrMsgTxt(errmsg);
    }
    }

    /* 7th Input. An array of at least length numcols containing the upper
       bound on each of the variables.
    */
    if ((nrhs > 6) && (mxGetM(UB_IN) != 0) && (mxGetN(UB_IN) != 0)) {
      if (!mxIsNumeric(UB_IN)
	  || (mxGetNumberOfDimensions(UB_IN) > 2)
	  || (mrowsc != mxGetM(UB_IN))
	  || (mxGetN(UB_IN) != 1)
	  || mxIsComplex(UB_IN)
	  || ((ub = mxGetPr(UB_IN)) == NULL)
	  ) {
	sprintf(errmsg,"UB must be a real valued %d by 1 column vector.",
		mrowsc);
	mexErrMsgTxt(errmsg);
      }
    }

    /* 8th Input. A column array containing the types of the variables.
     */
    if ((nrhs > 7) && (mxGetM(VARTYPE_IN) != 0)
	&& (mxGetN(VARTYPE_IN) != 0)) {
      if (!mxIsChar(VARTYPE_IN)
	  || (mxGetNumberOfDimensions(VARTYPE_IN) > 2)
	  || (mxGetM(VARTYPE_IN) != mrowsc)
	  || (mxGetN(VARTYPE_IN) != 1)
	  || mxIsComplex(VARTYPE_IN)
	  ) {
	sprintf(errmsg, "VARTYPE must be a char valued %d by 1 column vector.", mrowsc);
	mexErrMsgTxt(errmsg);
      } else {
	/*int i,size;*/