Implementation of subroutines

Solver library calls

The first two entries of the list of subroutines concern the handling of module parameters, with the exception of myxp_verbose that is a setting for the module itself, all other parameters are straightforward mappings of solver control parameters published by the solver library.

/**** Getting a control parameter ****/
static int slv_lc_getpar(XPRMcontext ctx,void *libctx)
{
 s_slvctx *slctx;
 int n;
 double r;

 slctx=libctx;
 n=XPRM_POP_INT(ctx);
 switch(n)
 {
  case 0:
   XPRM_PUSH_INT(ctx,(slctx->options&OPT_VERBOSE)?1:0);
   break;
  case 1:
   XPRSgetintcontrol(SLVCTX2PB(slctx)->xpb,XPRS_MAXTIME,&n);
   XPRM_PUSH_INT(ctx,n);
   break;
  case 2:
   XPRSgetintattrib(SLVCTX2PB(slctx)->xpb,XPRS_LPSTATUS,&n);
   XPRM_PUSH_INT(ctx,n);
   break;
  case 3:
   XPRSgetdblattrib(SLVCTX2PB(slctx)->xpb,XPRS_LPOBJVAL,&r);
   XPRM_PUSH_REAL(ctx,r);
   break;
  default:
   mm->dispmsg(ctx,"myxprs: Wrong control parameter number.\n");
   return XPRM_RT_ERROR;
 }
 return XPRM_RT_OK;
}

/**** Setting a control parameter ****/
static int slv_lc_setpar(XPRMcontext ctx,void *libctx)
{
 s_slvctx *slctx;
 int n;

 slctx=libctx;
 n=XPRM_POP_INT(ctx);
 switch(n)
 {
  case 0:
    slctx->options=XPRM_POP_INT(ctx)?(slctx->options|OPT_VERBOSE):(slctx->options&~OPT_VERBOSE);
    break;
  case 1:
    XPRSsetintcontrol(SLVCTX2PB(slctx)->xpb,XPRS_MAXTIME,XPRM_POP_INT(ctx));
    break;
  default:
    mm->dispmsg(ctx,"myxprs: Wrong control parameter number.\n");
    return XPRM_RT_ERROR;
 }

 return XPRM_RT_OK;
} 

The subroutine getprobstat exposes the Mosel problem status at the model level (this status value is populated after every solver run, see implementation of function slv_optim below).

static int slv_lc_getpstat(XPRMcontext ctx,void *libctx)
{
 XPRM_PUSH_INT(ctx,mm->getprobstat(ctx)&XPRM_PBRES);
 return RT_OK;
} 

The two module functions implementing the minimize and maximize subroutines map to the same function slv_optim.

static int slv_lc_maxim(XPRMcontext ctx,void *libctx)
{
 XPRMlinctr obj;

 obj=XPRM_POP_REF(ctx);
 return slv_optim(ctx,(s_slvctx *)libctx,OBJ_MAXIMIZE,obj);
}

static int slv_lc_minim(XPRMcontext ctx,void *libctx)
{
 XPRMlinctr obj;

 obj=XPRM_POP_REF(ctx);
 return slv_optim(ctx,(s_slvctx *)libctx,OBJ_MINIMIZE,obj);
}

The function slv_optim first clears any existing solution information, it then generates and loads the matrix representation of the problem into the solver and starts the actual solving process. After termination of the solver run it retrieves problem status information in order to populate Mosel's problem status flag.

static int slv_optim(XPRMcontext ctx, s_slvctx *slctx, int objsense, XPRMlinctr obj)
{
 int c,i;
 s_slvpb *slpb;
 int result;
 double objval;

 slpb=SLVCTX2PB(slctx);
 slpb->saved_ctx=ctx;	/* Save current context for callbacks */
 slv_clearsol(ctx,slpb);

 /* Call NI function 'loadmat' to generate and load the matrix */
 if(mm->loadmat(ctx,obj,NULL,MM_MAT_FORCE,&xpress,slpb)!=0)
 {
  mm->dispmsg(ctx,"myxprs: loadprob failed.\n");
  slpb->saved_ctx=NULL;
  return RT_ERROR;
 }

 /* Set optimization direction */
 XPRSchgobjsense(slpb->xpb,
  (objsense==OBJ_MINIMIZE)?XPRS_OBJ_MINIMIZE:XPRS_OBJ_MAXIMIZE);

 mm->setprobstat(ctx,XPRM_PBSOL,0);	/* Solution available for callbacks */
 if(!slpb->is_mip)
 {
  /* Solve an LP problem */
  c=XPRSlpoptimize(slpb->xpb,"");
  if(c!=0)
  {
   mm->dispmsg(ctx,"myxprs: optimisation failed.\n");
   slpb->saved_ctx=NULL;
   return RT_ERROR;
  }

  /* Retrieve solution status */
  XPRSgetintattrib(slpb->xpb,XPRS_PRESOLVESTATE,&i);
  if(i&128)
  {
   XPRSgetdblattrib(slpb->xpb,XPRS_LPOBJVAL,&objval);
   result=XPRM_PBSOL;
  }
  else
  {
   objval=0;
   result=0;
  }
  XPRSgetintattrib(slpb->xpb,XPRS_LPSTATUS,&i);
  switch(i)
  {
   case XPRS_LP_OPTIMAL:        result|=XPRM_PBOPT; break;
   case XPRS_LP_INFEAS:         result|=XPRM_PBINF; break;
   case XPRS_LP_CUTOFF:         result|=XPRM_PBOTH; break;
   case XPRS_LP_UNFINISHED:     result|=XPRM_PBUNF; break;
   case XPRS_LP_UNBOUNDED:      result|=XPRM_PBUNB; break;
   case XPRS_LP_CUTOFF_IN_DUAL: result|=XPRM_PBOTH; break;
   case XPRS_LP_UNSOLVED:       result|=XPRM_PBOTH; break;
  }
 }
 else
 {
  /* Solve an MIP problem */
  c=XPRSmipoptimize(slpb->xpb,"");
  if(c!=0)
  {
   mm->dispmsg(ctx,"myxprs: optimization failed.\n");
   slpb->saved_ctx=NULL;
   return RT_ERROR;
  }

  /* Retrieve solution status */
  XPRSgetintattrib(slpb->xpb,XPRS_MIPSTATUS,&i);
  switch(i)
  {
   case XPRS_MIP_LP_NOT_OPTIMAL:
      objval=0;
      result=XPRM_PBUNF;
      break;
   case XPRS_MIP_LP_OPTIMAL:
      objval=0;
      result=XPRM_PBUNF;
      break;
   case XPRS_MIP_NO_SOL_FOUND:	/* Search incomplete: no solution */
      objval=0;
      result=XPRM_PBUNF;
      break;
   case XPRS_MIP_SOLUTION:	/* Search incomplete: there is a solution */
      XPRSgetdblattrib(slpb->xpb,XPRS_MIPOBJVAL,&objval);
      result=XPRM_PBUNF|XPRM_PBSOL;
      slpb->have|=HAVEMIPSOL;
      break;
   case XPRS_MIP_INFEAS:	/* Search complete: no solution */
      objval=0;
      result=XPRM_PBINF;
      break;
   case XPRS_MIP_OPTIMAL:	/* Search complete: best solution available */
      XPRSgetdblattrib(slpb->xpb,XPRS_MIPOBJVAL,&objval);
      result=XPRM_PBSOL|XPRM_PBOPT;
      slpb->have|=HAVEMIPSOL;
      break;
   case XPRS_MIP_UNBOUNDED:
      objval=0;
      result=XPRM_PBUNB;
      break;
  }
  if(!(result&XPRM_PBSOL))
  {
   /* If no MIP solution try to get an LP solution */
   XPRSgetintattrib(slpb->xpb,XPRS_PRESOLVESTATE,&i);
   if(i&128)
   {
    XPRSgetdblattrib(slpb->xpb,XPRS_LPOBJVAL,&objval);
    result|=XPRM_PBSOL;
   }
  }
 }

 /* Record solution status and objective value */
 mm->setprobstat(ctx,result,objval);
 slpb->saved_ctx=NULL;
 return 0;
} 

Implementation of MIP solver interface functions

The following functions implement the MIP solver interface functions ('loadmat', 'clearsol', 'getsol_v', 'getsol_c') that are communicated to the NI routine loadmat via the MIP solution interface structure xpress (for its definition see Section Module context).

The 'loadmat' function loads a matrix that is held in Mosel structures into the LP or MIP solver.

static int slv_loadmat(XPRMcontext ctx, void *mipctx, mm_matrix *m)
{
 s_slvpb *slpb;
 s_slvctx *slctx;
 char pbname[80];
 int c,r;

 slpb=mipctx;
 slctx=slpb->slctx;
 slv_clearsol(ctx,slpb);
 slpb->is_mip=(m->ngents>0) || (m->nsos>0);

 sprintf(pbname,"xpb%p",slpb);
 if(slpb->is_mip)
  r=XPRSloadglobal(slpb->xpb,pbname,m->ncol,m->nrow,
		m->qrtype,m->rhs,m->range,m->obj,
		m->mstart,NULL,m->mrwind,m->dmatval,m->dlb,m->dub,
		m->ngents,m->nsos,m->qgtype,m->mgcols,m->mplim,m->qstype,
		m->msstart,m->mscols,m->dref);
 else
  r=XPRSloadlp(slpb->xpb,pbname,m->ncol,m->nrow,
		m->qrtype,m->rhs,m->range,m->obj,
		m->mstart,NULL,m->mrwind,m->dmatval,m->dlb,m->dub);

 /* Objective constant term */
 if(!r) { c=-1;r=XPRSchgobj(slpb->xpb,1, &c, &(m->fixobj)); }

 return r;
} 

The 'clearsol' function frees up solution information held in the solver problem interface structures.

static void slv_clearsol(XPRMcontext ctx, void *mipctx)
{
 s_slvpb *slpb;

 slpb=mipctx;
 Free(&(slpb->solval));
 Free(&(slpb->dualval));
 Free(&(slpb->rcostval));
 Free(&(slpb->slackval));
 slpb->have=0;
}

/**** Free + reset memory ****/
static void Free(void *ad)
{
 free(*(void **)ad);
 *(void **)ad=NULL;
}

The 'getsol_v' and 'getsol_c' routines serve to retrieve the solution values for decision variables and constraints from the solver. These implementations retrieve the entire arrays at once and any subsequent calls return the information saved in the solver interface structures.

/**** Solution information for decision variables ****/
static double slv_getsol_v(XPRMcontext ctx, void *mipctx, int what, int col)
{
 s_slvpb *slpb;

 slpb=mipctx;
 if(what)
 {
  if(!(slpb->have&HAVERCS))
  {
   if(slpb->rcostval==NULL)
   {
    int ncol;

    XPRSgetintattrib(slpb->xpb,XPRS_ORIGINALCOLS,&ncol);
    if((slpb->rcostval=malloc(ncol*sizeof(double)))==NULL)
    {
     mm->dispmsg(ctx,"myxprs: Out of memory error.\n");
     return 0;
    }
   }
   if(slpb->have&HAVEMIPSOL)	/* No rcost for a MIP => 0 */
   {
    int ncol;

    XPRSgetintattrib(slpb->xpb,XPRS_ORIGINALCOLS,&ncol);
    memset(slpb->rcostval,0,ncol*sizeof(double));
   }
   else
    XPRSgetlpsol(slpb->xpb,NULL,NULL,NULL,slpb->rcostval);
   slpb->have|=HAVERCS;
  }
  return slpb->rcostval[col];
 }
 else
 {
  if(!(slpb->have&HAVESOL))
  {
   if(slpb->solval==NULL)
   {
    int ncol;

    XPRSgetintattrib(slpb->xpb,XPRS_ORIGINALCOLS,&ncol);
    if((slpb->solval=malloc(ncol*sizeof(double)))==NULL)
    {
     mm->dispmsg(ctx,"myxprs: Out of memory error.\n");
     return 0;
    }
   }
   if(slpb->have&HAVEMIPSOL)
    XPRSgetmipsol(slpb->xpb,slpb->solval,NULL);
   else
    XPRSgetlpsol(slpb->xpb,slpb->solval,NULL,NULL,NULL);
   slpb->have|=HAVESOL;
  }
  return slpb->solval[col];
 }
} 

/**** Solution information for linear constraints ****/
static double slv_getsol_c(XPRMcontext ctx, void *mipctx, int what, int row)
{
 s_slvpb *slpb;

 slpb=mipctx;
 if(what)
 {
  if(!(slpb->have&HAVEDUA))
  {
   if(slpb->dualval==NULL)
   {
    int nrow;

    XPRSgetintattrib(slpb->xpb,XPRS_ORIGINALROWS,&nrow);
    if((slpb->dualval=malloc(nrow*sizeof(double)))==NULL)
    {
     mm->dispmsg(ctx,"myxprs: Out of memory error.\n");
     return 0;
    }
   }
   if(slpb->have&HAVEMIPSOL)	/* No dual for a MIP => 0 */
   {
    int nrow;

    XPRSgetintattrib(slpb->xpb,XPRS_ORIGINALROWS,&nrow);
    memset(slpb->dualval,0,nrow*sizeof(double));
   }
   else
    XPRSgetlpsol(slpb->xpb,NULL,NULL,slpb->dualval,NULL);
   slpb->have|=HAVEDUA;
  }
  return slpb->dualval[row];
 }
 else
 {
  if(!(slpb->have&HAVESLK))
  {
   if(slpb->slackval==NULL)
   {
    int nrow;

    XPRSgetintattrib(slpb->xpb,XPRS_ORIGINALROWS,&nrow);
    if((slpb->slackval=malloc(nrow*sizeof(double)))==NULL)
    {
     mm->dispmsg(ctx,"myxprs: Out of memory error.\n");
     return 0;
    }
   }
   if(slpb->have&HAVEMIPSOL)
    XPRSgetmipsol(slpb->xpb,NULL,slpb->slackval);
   else
    XPRSgetlpsol(slpb->xpb,NULL,slpb->slackval,NULL,NULL);
   slpb->have|=HAVESLK;
  }
  return slpb->slackval[row];
 }
} 

Implementation of services

The RESET service function creates a new module context if none is provided in the argument, on a subsequent call where the module context argument is populated this context will be released after freeing all data structures that may have been created via the solver library.

/**** Reset the myxprs interface for a run ****/
static void *slv_reset(XPRMcontext ctx, void *libctx)
{
 s_slvctx *slctx;

 /* End of execution: release context */
 if(libctx!=NULL)
 {
  slctx=libctx;

  /* Release all remaining problems */
  while(slctx->probs!=NULL)
  {
   slv_pb_delete(ctx,slctx,slctx->probs,-1);
  }
  free(slctx);
  return NULL;
 }
 else
 {
  /* Begin of execution: create context */
  if((slctx=malloc(sizeof(s_slvctx)))==NULL)
  {
   mm->dispmsg(ctx,"myxprs: Out of memory error.\n");
   return NULL;
  }
  memset(slctx,0,sizeof(s_slvctx));

  /* Record the problem ID of our problem type */
  mm->gettypeprop(ctx,mm->findtypecode(ctx,"mpproblem.mxp"),XPRM_TPROP_PBID,
  						(XPRMalltypes*)&(slctx->pbid));
  return (void *)slctx;
 }
}

The UNLOAD service terminates the solver library (frees up the licence) that has been initialized from the module initialization.

/**** Called when unloading the library ****/
static void slv_quitlib(void)
{
 if(mm!=NULL)
 {
  XPRSfree();
  mm=NULL;
 }
} 

The implementation of the module parameter access services PARAM and PARLST is similar to what we have seen for other modules (e.g. see Section Services related to parameters).

/**** Find a control parameter ****/
static int slv_findparam(const char *name,int *type,int why,XPRMcontext ctx,
								void *libctx)
{
 int n;

 for(n=0;n<SLV_NBPARAM;n++)
 {
  if(strcmp(name,myxprsparams[n].name)==0)
  {
   *type=myxprsparams[n].type;
   return n;
  }
 }
 return -1;
}

/**** Return the next parameter for enumeration ****/
static void *slv_nextparam(void *ref,const char **name,const char **desc,
								int *type)
{
 long cst;

 cst=(long)ref;
 if((cst<0)||(cst>=SLV_NBPARAM))
  return NULL;
 else
 {
  *name=myxprsparams[cst].name;
  *type=myxprsparams[cst].type;
  *desc=NULL;
  return (void *)(cst+1);
 }
} 

Handling optimization problems

Each Mosel model creates a default optimization problem of type mpproblem holding the constraints that are defined in the model. Further (sub)problems can be defined explicitly by the model developer, such as in the example shown at the beginning of this chapter (Section Example).

A solver module needs to implement an extension to the mpproblem type. Typically, the underlying data structure will include a reference to the solver problem representation, some status flags and structures to store solution information (see definition in Section List of types above). For our Xpress Optimizer example we have implemented the type handling routines to create, delete, and copy optimization problems. If the underlying solver can only handle a single problem, the implementation of the 'create' routine should prevent the creation of more than one problem and a 'copy' routine is most likely not required.

The following implementation of a problem creation routine for Xpress Optimizer creates the Optimizer problem, redirects the Optimizer output onto Mosel and it also defines some logging callbacks in order to intercept a program interruption.

/**** Create a new "problem" ****/
static void *slv_pb_create(XPRMcontext ctx, void *libctx, void *toref, int type)
{
 s_slvctx *slctx;
 s_slvpb *slpb;
 int i;

 slctx=libctx;
 if((slpb=malloc(sizeof(s_slvpb)))==NULL)
 {
  mm->dispmsg(ctx,"myxprs: Out of memory error.\n");
  return NULL;
 }
 memset(slpb,0,sizeof(s_slvpb));
 i=XPRScreateprob(&(slpb->xpb));
 if((i!=0)&&(i!=32))
 {
  mm->dispmsg(ctx,"myxprs: I cannot create the problem.\n");
  free(slpb);
  return NULL;
 }
 slpb->slctx=slctx;
 /* Redirect solver messages to the Mosel streams */
 XPRSaddcbmessage(slpb->xpb,slvcb_output,slpb,0);
 XPRSsetintcontrol(slpb->xpb,XPRS_OUTPUTLOG,1);

 /* Define log callbacks to report program interruption */
 XPRSaddcblplog(slpb->xpb,(void*)slvcb_stopxprs,slpb,0);
 XPRSaddcbcutlog(slpb->xpb,(void*)slvcb_stopxprs,slpb,0);
 XPRSaddcbgloballog(slpb->xpb,(void*)slvcb_stopxprs,slpb,0);
 XPRSaddcbbarlog(slpb->xpb,(void*)slvcb_stopxprs,slpb,0);

 if(slctx->probs!=NULL)
 {
  slpb->next=slctx->probs;
  slctx->probs->prev=slpb;
 }
 /* else we are creating the master problem */

 slctx->probs=slpb;
 return slpb;
} 

The problem deletion routine needs to update the list of problems saved in the solver problem interface structure.

/**** Delete a "problem" ****/
static void slv_pb_delete(XPRMcontext ctx,void *libctx,void *todel,int type)
{
 s_slvctx *slctx;
 s_slvpb *slpb;

 slctx=libctx;
 slpb=todel;
 slv_clearsol(ctx,slpb);
 XPRSdestroyprob(slpb->xpb);
 if(slpb->next!=NULL)	/* Last in list */
  slpb->next->prev=slpb->prev;
 if(slpb->prev==NULL)	/* First in list */
  slctx->probs=slpb->next;
 else
  slpb->prev->next=slpb->next;
 free(slpb);
} 

A problem is copied without duplicating the solution information.

/**** Copy/reset/append problems: simply clear data of the destination ****/
static int slv_pb_copy(XPRMcontext ctx,void *libctx,void *toinit,void *src,int ust)
{
 s_slvpb *slpb;

 slpb=toinit;
 if(XPRM_CPY(ust)<XPRM_CPY_APPEND) slv_clearsol(ctx,slpb);
 return 0;
}