/*---------------------------------------------------------------------- * * Copyright (C) 2001-2005, Nicolo' Giorgetti, All rights reserved. * E-mail: <giorgetti@dii.unisi.it>. * * 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-14void 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;*/