glpkmex.c

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

	size = mxGetNumberOfElements(VARTYPE_IN)+1;
	/*   Allocate enough memory to hold the converted string.
	 */
	vartype   = mxCalloc(size, sizeof (char));
	vartype2 = mxCalloc(size, sizeof (int));

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

	/* checking if the input is made only of C I */
	for (i = 0; i < size-1 ; i++){
	  if ((vartype[i] != 'C') && (vartype[i] != 'I'))
	    mexErrMsgTxt("VARTYPE must contain only 'C' or 'I'");
	  if(vartype[i]=='I'){
	    isMIP=1;
	    vartype2[i]=LPX_IV;
	  }else{
	    vartype2[i]=LPX_CV;
	  }
	}
      }
    } else mexWarnMsgTxt ("Omitting VARTYPE you are assuming all variables are continuous");

    /*
           9th Input. A structure containing the control parameters.
	*/
	if ((nrhs > 8) && ((nfields=mxGetNumberOfFields(PARAM)) !=0)) {
	  int *idtmp=NULL;
	  if(!mxIsStruct(PARAM)){
	    mexErrMsgTxt("PARAM must be a structure !");
	  }


	  /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
	  /* Integer parameters                                             */
	  /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

	  /* level of messages output by the solver */
	  if((mxtmp=mxGetField(PARAM,0,"msglev")) != NULL){
	    rdtmp=mxGetPr(mxtmp);
	    if((*rdtmp != 0) && (*rdtmp != 1) && (*rdtmp != 2) && (*rdtmp != 3)){
	      sprintf(errmsg,"'msglev' parameter must be only:\n\t0 - no output,\n\t1 - error messages only),\n\t2 - normal output,\n\t3 - full output [default]");
	      mexErrMsgTxt(errmsg);
	    }
	    lpxIntParam[0]=(int) *rdtmp;
	  }

	    /* scaling option */
	    if((mxtmp=mxGetField(PARAM,0,"scale")) != NULL){
	      rdtmp=mxGetPr(mxtmp);
	      if((*rdtmp != 0) && (*rdtmp != 1) && (*rdtmp != 2)){
				sprintf(errmsg,"'scale' parameter must be only:\n\t0 - no scaling,\n\t1 - equilibration scaling,\n\t2 - geometric mean scaling");
				mexErrMsgTxt(errmsg);
	      }
	      lpxIntParam[1]=(int) *rdtmp;
	    }

	    /* Dual dimplex option */
	    if((mxtmp=mxGetField(PARAM,0,"dual")) != NULL){
	      rdtmp=mxGetPr(mxtmp);
	      if((*rdtmp != 0) && (*rdtmp != 1)){
				sprintf(errmsg,"'dual' parameter must be only:\n\t0 - do not use the dual simplex [default],\n\t1 - use dual simplex");
				mexErrMsgTxt(errmsg);
	      }
	      lpxIntParam[2]=(int) *rdtmp;
	    }

	    /* pricing option */
	    if((mxtmp=mxGetField(PARAM,0,"price"))  != NULL){
	      rdtmp=mxGetPr(mxtmp);
	      if((*rdtmp != 0) && (*rdtmp != 1)){
				sprintf(errmsg,"'price' parameter must be only:\n\t0 - textbook pricing,\n\t1 - steepest edge pricing [default]");
				mexErrMsgTxt(errmsg);
	      }
	      lpxIntParam[3]=(int) *rdtmp;
	    }

	    /* solution rounding option */
	    if((mxtmp=mxGetField(PARAM,0,"round")) != NULL){
	      rdtmp=mxGetPr(mxtmp);
	      if((*rdtmp != 0) && (*rdtmp != 1)){
				sprintf(errmsg,"'round' parameter must be only:\n\t0 - report all primal and dual values [default],\n\t1 - replace tiny primal and dual values by exact zero");
				mexErrMsgTxt(errmsg);
	      }
	      lpxIntParam[4]=(int) *rdtmp;
	    }

	    /* simplex iterations limit */
	    if((mxtmp=mxGetField(PARAM,0,"itlim")) != NULL){
	      rdtmp=mxGetPr(mxtmp); lpxIntParam[5]=(int) *rdtmp; }

	    /* Simplex iterations count */
	    if((mxtmp=mxGetField(PARAM,0,"itcnt")) != NULL){
	      rdtmp=mxGetPr(mxtmp); lpxIntParam[6]=(int) *rdtmp; }

	    /* Output frequency, in iterations */
	    if((mxtmp=mxGetField(PARAM,0,"outfrq")) != NULL){
	      rdtmp=mxGetPr(mxtmp); lpxIntParam[7]=(int) *rdtmp; }

	    /* Branching heuristic option  */
	    if((mxtmp=mxGetField(PARAM,0,"branch")) != NULL){
	      rdtmp=mxGetPr(mxtmp);
	      if((*rdtmp != 0) && (*rdtmp != 1) && (*rdtmp != 2)){
				sprintf(errmsg,"'branch' parameter must be only (for MIP only):\n\t0 - branch on the first variable,\n\t1 - branch on the last variable,\n\t2 - branch using a heuristic by Driebeck and Tomlin [default]");
				mexErrMsgTxt(errmsg);
	      }
	      lpxIntParam[14]=(int) *rdtmp;
	    }

	    /* Backtracking heuristic option  */
	    if((mxtmp=mxGetField(PARAM,0,"btrack")) != NULL){
	      rdtmp=mxGetPr(mxtmp);
	      if((*rdtmp != 0) && (*rdtmp != 1) && (*rdtmp != 2)){
				sprintf(errmsg,"'btrack' parameter must be only (for MIP only):\n\t0 - depth first search,\n\t1 - breadth first search,\n\t2 - backtrack using the best projection heuristic");
				mexErrMsgTxt(errmsg);
	      }
	      lpxIntParam[15]=(int) *rdtmp;
	    }

      if((mxtmp=mxGetField(PARAM,0,"presol")) != NULL){
         rdtmp=mxGetPr(mxtmp);
         if((*rdtmp != 0) && (*rdtmp != 1)){
            sprintf(errmsg,"'presol' parameter must be only:\n\t0 - LP presolver is ***NOT*** used,\n\t1 - LP presol is used");
            mexErrMsgTxt(errmsg);
         }
         lpxIntParam[16]=(int) *rdtmp;
      }

      /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
		/* Real parameters                                                */
		/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

		/* ratio test option */
		if((mxtmp=mxGetField(PARAM,0,"relax")) != NULL){
		   rdtmp=mxGetPr(mxtmp); lpxRealParam[0]=*rdtmp; }

		/* relative tolerance used to check if the current basic solution
		   is primal feasible */
		if((mxtmp=mxGetField(PARAM,0,"tolbnd")) != NULL){
		   rdtmp=mxGetPr(mxtmp); lpxRealParam[1]=*rdtmp; }

		/* absolute tolerance used to check if the current basic solution
			is dual feasible */
		if((mxtmp=mxGetField(PARAM,0,"toldj")) != NULL){
			rdtmp=mxGetPr(mxtmp); lpxRealParam[2]=*rdtmp; }

		/* relative tolerance used to choose eligible pivotal elements of
			the simplex table in the ratio test */
		if((mxtmp=mxGetField(PARAM,0,"tolpiv")) != NULL){
			rdtmp=mxGetPr(mxtmp); lpxRealParam[3]=*rdtmp; }

		if((mxtmp=mxGetField(PARAM,0,"objll")) != NULL){
			rdtmp=mxGetPr(mxtmp); lpxRealParam[4]=*rdtmp; }

		if((mxtmp=mxGetField(PARAM,0,"objul")) != NULL){
		  	rdtmp=mxGetPr(mxtmp); lpxRealParam[5]=*rdtmp; }

		if((mxtmp=mxGetField(PARAM,0,"tmlim")) != NULL){
		  rdtmp=mxGetPr(mxtmp); lpxRealParam[6]=*rdtmp; }

		if((mxtmp=mxGetField(PARAM,0,"outdly")) != NULL){
		  rdtmp=mxGetPr(mxtmp); lpxRealParam[7]=*rdtmp; }

		if((mxtmp=mxGetField(PARAM,0,"tolint")) != NULL){
		  rdtmp=mxGetPr(mxtmp); lpxRealParam[8]=*rdtmp; }

		if((mxtmp=mxGetField(PARAM,0,"tolobj")) != NULL){
		  rdtmp=mxGetPr(mxtmp); lpxRealParam[9]=*rdtmp; }

	}

	/* 10th Input. If the problem is a LP problem you may select which solver
	   use: RSM (Revised Simplex Method) or IPM (Interior Point Method).
	   If the problem is a MIP problem this field will be ignored.
	*/
	if ((nrhs > 9) && (mxGetM(SOLVER_IN) != 0) && (mxGetN(SOLVER_IN) != 0)) {
	  if (!mxIsNumeric(SOLVER_IN)
	      || (mxGetNumberOfDimensions(SOLVER_IN) > 2)
	      || (1 != mxGetM(SOLVER_IN))
	      || (mxGetN(SOLVER_IN) != 1)
	      || mxIsComplex(SOLVER_IN)
	      || ((tmp = mxGetPr(SOLVER_IN)) == NULL)
	      ) {
	    mexErrMsgTxt("SOLVER_IN must be a real valued scalar.");
	  } else {
	    lpsolver = (int) tmp[0];
	  }
	}

	/* 11th Input. Saves a copy of the problem if SAVE<>0.
	 */
	if ((nrhs > 10) && (mxGetM(SAVE_IN) != 0) && (mxGetN(SAVE_IN) != 0)) {
	  if (!mxIsNumeric(SAVE_IN)
	      || (mxGetNumberOfDimensions(SAVE_IN) > 2)
	      || (1 != mxGetM(SAVE_IN))
	      || (mxGetN(SAVE_IN) != 1)
	      || mxIsComplex(SAVE_IN)
	      || ((tmp = mxGetPr(SAVE_IN)) == NULL)
	      ) {
	    mexErrMsgTxt("SAVE must be a real valued scalar.");
	  } else {
	    save_pb = (tmp[0] != 0);
	  }
	}


	if (nrhs > 11) mexWarnMsgTxt("Extra parameters ignored.");


	/* ---- Set default values ---- */
	/* CTYPE argument, default: upper bound */
	if (ctype == NULL) {
	  /*int i;*/
	  ctype = mxCalloc(mrowsA, sizeof (char));

	  for (i = 0; i < mrowsA; i++) ctype[i] = 'U';
	}
	
	/* B argument: remove possible infinity values */
	bcopy=(double *)mxCalloc(mrowsA,sizeof(double));
	for(i=0; i< mrowsA; i++){
		if (b[i]==-mxGetInf()) bcopy[i]=-1e12;
		if (b[i]==mxGetInf()) bcopy[i]=1e12;
		else bcopy[i]=b[i];
	}	
	
	/*LB argument, default: Free */
	freeLB=(int *) mxCalloc(mrowsc, sizeof(int));
	freeUB=(int *) mxCalloc(mrowsc, sizeof(int));
	if (lb == NULL) {
	  flagLB=1;
	  /*int i;*/
	  lb = mxCalloc(mrowsc, sizeof (double));
	  for (i = 0; i < mrowsc; i++){
	  	 lb[i] = -mxGetInf();
	  	 freeLB[i]=1;
	  }
	}else{
	  /*int i;*/
		for(i=0;i<mrowsc;i++){
			if(lb[i]==-mxGetInf()) freeLB[i]=1;
			else freeLB[i]=0;
		}
	}
	/*UB argument, default: Free */
	if (ub == NULL) {
	  flagUB=1;
	  /*int i;*/
	  ub = mxCalloc(mrowsc, sizeof (double));
	  for (i = 0; i < mrowsc; i++){
	  	ub[i] = mxGetInf();
	  	freeUB[i]=1;
	  }
	}else{
	  /*int i;*/
		for(i=0;i<mrowsc;i++){
			if(ub[i]==mxGetInf()) freeUB[i]=1;
			else freeUB[i]=0;
		}
	}
	/*VARTYPE argument, default: continuous */
	if (vartype == NULL) {
	  /*int i;*/
	  vartype2 = mxCalloc(mrowsc, sizeof (int));

	  for (i = 0; i < mrowsc; i++) vartype2[i] = LPX_CV;
	}

	extranames=mxCalloc(4,sizeof(*extranames));
	extranames[0]="lambda";
	extranames[1]="redcosts";
	extranames[2]="time";
	extranames[3]="memory";
	/* Create a matrices for the return arguments
	 */
	XMIN_OUT   = mxCreateDoubleMatrix(mrowsc, 1, mxREAL);
	FMIN_OUT   = mxCreateDoubleMatrix(1, 1, mxREAL);
	STATUS_OUT = mxCreateDoubleMatrix(1, 1, mxREAL);
	
	EXTRA_OUT  = mxCreateStructMatrix(1, 1, 4, extranames);
	mxlambda   = mxCreateDoubleMatrix(mrowsA, 1, mxREAL);
	mxredcosts  = mxCreateDoubleMatrix(mrowsc, 1, mxREAL);
	mxtime     = mxCreateDoubleMatrix(1, 1, mxREAL);
	mxmem      = mxCreateDoubleMatrix(1, 1, mxREAL);
	

	/* Assign pointers to the output parameters
	 */
	xmin   = mxGetPr(XMIN_OUT);
	fmin   = mxGetPr(FMIN_OUT);
	status = mxGetPr(STATUS_OUT);
	lambda = mxGetPr(mxlambda);
	redcosts= mxGetPr(mxredcosts);
	time   = mxGetPr(mxtime);
	mem    = mxGetPr(mxmem);

	jmpret = setjmp( mark );
   if (jmpret==0){
	error=glpk(sense,mrowsc,mrowsA,c,nz,rn,cn,a,bcopy,ctype,
		    freeLB,lb,freeUB,ub,vartype2,isMIP,lpsolver,save_pb,
		    xmin,fmin,status,lambda,redcosts, time,mem);
	}

	mxSetField(EXTRA_OUT,0,extranames[0],mxlambda);
	mxSetField(EXTRA_OUT,0,extranames[1],mxredcosts);
	mxSetField(EXTRA_OUT,0,extranames[2],mxtime);
	mxSetField(EXTRA_OUT,0,extranames[3],mxmem);

	if(ctype != NULL) mxFree(ctype);
	if(vartype != NULL) mxFree(vartype);
	if(vartype2 != NULL) mxFree(vartype2);
	if(flagLB) mxFree(lb);
	if(flagUB) mxFree(ub);
	mxFree(extranames);
	mxFree(rn);
	mxFree(cn);
	mxFree(a);
	mxFree(freeLB);
	mxFree(freeUB);
	mxFree(bcopy);
	return;
}