Example: implementing an interface to a QPQC solver

The example implementation described here can be considered as an extension with handling of nonlinearities of the example documented in the chapter Implementing an LP/MIP solver interface of the Mosel NI User Guide. For solver-specific features that are not directly related to the handling of nonlinear constraints—such as the implementation of callback functions—the reader is therefore referred to this document.

Mosel model using myqxprs

The following example shows how to calculate the shape of a hanging chain with N chainlinks and fixed endpoints. In this example the objective function is linear and the constraints contain nonlinear (quadratic) terms. The modelling features that are provided through the example solver module myqxprs are highlighted. Features such as the setting of initial values for decision variables, the retrieval of solution values, or the type nlctr for nonlinear constraints are provided by the module mmnl that is included by our solver module (see Chapter mmnl of the Mosel Language Reference Manual for the documentation of mmnl).

model "catenary - myqxprs version"
 uses "myqxprs"

 parameters
   N = 100                       ! Number of chainlinks
   L = 1                         ! Difference in x-coordinates of endlinks
   H = 2*L/N                     ! Length of each link
 end-parameters

 declarations
   RN = 0..N
   x: array(RN) of mpvar         ! x-coordinates of endpoints of chainlinks
   y: array(RN) of mpvar         ! y-coordinates of endpoints of chainlinks
   PotentialEnergy: linctr       ! Objective function
   Link: array(range) of nlctr   ! Constraints
 end-declarations

 forall(i in RN) x(i) is_free
 forall(i in RN) y(i) is_free

! Objective: minimise the potential energy
 PotentialEnergy:= sum(j in 1..N) (y(j-1)+y(j))/2

! Bounds: positions of endpoints
 x(0) = 0; y(0) = 0; x(N) = L; y(N) = 0

! Constraints: positions of chainlinks
 forall(j in 1..N)
   Link(j):= (x(j)-x(j-1))^2+(y(j)-y(j-1))^2 <= H^2

! Setting start values
 forall(j in RN) setinitval(x(j), j*L/N)
 forall(j in RN) setinitval(y(j), 0)

! Configuration of the solver
 setparam("myxp_verbose", true)

! Solve the problem
 minimise(PotentialEnergy)

! Solution reporting
 if getprobstat not in {MYXP_OPT, MYXP_UNF} then
   writeln("No solution available. Solver status: ", getprobstat)
 else
   writeln("Solution: ", getobjval)
   forall(j in RN)
     writeln(strfmt(getsol(x(j)),10,5), " ", strfmt(getsol(y(j)),10,5))
 end-if

end-model

Structures for passing information

A minimal implementation of a solver module based on mmnl needs to do the following:

• Modeling functionality:
  • define subroutines to start an optimization run and retrieve solution values
  • provide access to solver parameters
  • if supported by the solver, provide support for handling multiple problems
• NI functionality:
  • implement a reset and an unload service
  • implement module dependency services for mmnl
  • initialize the module and the required interface structures

The next few sections (List of subroutines – Module context) describe the structures that are required for exchanging information between Mosel and the external program (the solver library): these structures are communicated during the module initialization (see Initialization function). This is followed by an explanation of the implementation of the main module routines starting in Section Implementation of subroutines.

Beyond the minimal functionality of a solver module there are many other tasks that one might wish to expose at the Mosel language level, one such example is given in Section Generating names for matrix entries and matrix output with the writeprob matrix output routine and the related Mosel NI functionality for generating names based on model entity names. Another example would be the implementation of solver callback functions that allow users to interact with the problem via the Mosel language at specific points while the solver is running—please see Section Implementing a solver callback of the Mosel NI User Guide for the description of an implementation example.

List of subroutines

The minimal set of entries for the list of subroutines would be just the calls to minimization/maximization. Our implementation adds a function to retrieve the problem status information, alternative spelling for the optimization routines, and it also provides the access routines for module control parameters that are required for the implementation of solver parameters. The optimization routines have two different implementations in our example, one for linear expressions(of type linctr) and a second for nonlinear expressions (the type nlctr defined by mmnl).

static XPRMdsofct tabfct[]=
    {
     {"",XPRM_FCT_GETPAR,XPRM_TYP_NOT,0,NULL,slvlc_getpar},
     {"",XPRM_FCT_SETPAR,XPRM_TYP_NOT,0,NULL,slvlc_setpar},
     {"getprobstat",2000,XPRM_TYP_INT,0,NULL,slvlc_getpstat},
     {"minimise",2100,XPRM_TYP_NOT,1,"c",slvlc_minim},
     {"minimize",2100,XPRM_TYP_NOT,1,"c",slvlc_minim},
     {"maximise",2101,XPRM_TYP_NOT,1,"c",slvlc_maxim},
     {"maximize",2101,XPRM_TYP_NOT,1,"c",slvlc_maxim}
     {"minimise", 2120, XPRM_TYP_NOT,1,"|nlctr|", slv_lc_nminim},
     {"minimize", 2120, XPRM_TYP_NOT,1,"|nlctr|", slv_lc_nminim},
     {"maximise", 2121, XPRM_TYP_NOT,1,"|nlctr|", slv_lc_nmaxim},
     {"maximize", 2121, XPRM_TYP_NOT,1,"|nlctr|", slv_lc_nmaxim}
    }; 

List of parameters

In terms of an example, we provide access to a few controls of Xpress Optimizer, and the module also shows how to implement a verbosity flag, resulting in the following list of module parameters:

static struct                   /* Parameters published by this module */
    {
     char *name;
     int type;
    } myxprsparams[]=
    {
     {"myxp_verbose", XPRM_TYP_BOOL|XPRM_CPAR_READ|XPRM_CPAR_WRITE},
     {"myxp_maxtime", XPRM_TYP_INT|XPRM_CPAR_READ|XPRM_CPAR_WRITE},
     {"myxp_lpstatus", XPRM_TYP_INT|XPRM_CPAR_READ},
     {"myxp_lpobjval", XPRM_TYP_REAL|XPRM_CPAR_READ},
    }; 

The problem and LP status parameters return values are best implemented via module constants, such as:

static XPRMdsoconst tabconst[]=
    {
     XPRM_CST_INT("MYXP_INF",XPRM_PBINF),             /* Mosel status codes */
     XPRM_CST_INT("MYXP_OPT",XPRM_PBOPT),
     XPRM_CST_INT("MYXP_OTH",XPRM_PBOTH),
     XPRM_CST_INT("MYXP_UNF",XPRM_PBUNF),
     XPRM_CST_INT("MYXP_UNB",XPRM_PBUNB),
     XPRM_CST_INT("MYXP_LP_OPTIMAL",XPRS_LP_OPTIMAL), /* Solver status codes */
     XPRM_CST_INT("MYXP_LP_INFEAS",XPRS_LP_INFEAS),
     XPRM_CST_INT("MYXP_LP_CUTOFF",XPRS_LP_CUTOFF)
    }; 

List of types

The list of types has a single entry: a solver module needs to extend the Mosel type mpproblem with its own implementation.

static XPRMdsotyp tabtyp[]=
    {
     {"mpproblem.mqxp", 1, XPRM_DTYP_PROB|XPRM_DTYP_APPND, slv_pb_create,
      slv_pb_delete, NULL, NULL, slv_pb_copy}
    }; 

The following structure implements the problem type for our module, Mosel will maintain one instance of this type for each mpproblem object.

typedef struct SlvPb
    {
     struct SlvCtx *slctx;      /* Solver module context */
     XPRSprob xpb;
     int have;
     double *solval;            /* Structures for storing solution values */
     double *dualval;
     double *rcostval;
     double *slackval;
     XPRMcontext saved_ctx;     /* Mosel context (used by callbacks) */
     struct SlvPb *prev,*next;
    } s_slvpb;

A solver module context definition is shown below in Section Module context.

List of services

The services PARAM and PARLST are required for the handling of module parameters, RESET and UNLOAD manage the access to the solver library. Furthermore, we need to specify the dependency on functionality provided via the NI interface of the module mmnl via the DEPLST and IMPLST services.

static XPRMdsoserv tabserv[]=
    {
     {XPRM_SRV_PARAM, (void *)slv_findparam},
     {XPRM_SRV_PARLST, (void *)slv_nextparam},
     {XPRM_SRV_RESET, (void *)slv_reset},
     {XPRM_SRV_UNLOAD, (void *)slv_quitlib}
     {XPRM_SRV_DEPLST, (void *)slv_deplist},
     {XPRM_SRV_IMPLST, (void *)slv_implst}
    }; 

The first four entries in this list are implemented via module (C library) functions (see Section Implementation of services below for their immplementation), the module dependencies are stated via these structures:

static const char *slv_deplist[]={"mmnl",NULL};
static const char *slv_implst[]={"mmnl",NULL};

Module context

The module context holds the type ID for the extended mpproblem type, module options, and a list of references to the problems that have been created by this module:

typedef struct SlvCtx           /* A context for this module */
    {
     int pbid;                  /* ID of type "mpproblem.mqxp" */
     int options;               /* Runtime options */
     s_slvpb *probs;            /* List of created problems */
     void **nlctx;              /* Execution context of 'mmnl' */
    } s_slvctx; 

Interface structure

The interface structure that is passed to Mosel holds the definition of the four lists (constants, subroutines, types, and services).

static XPRMdsointer dsointer=
    {
     sizeof(tabconst)/sizeof(XPRMdsoconst), tabconst,
     sizeof(tabfct)/sizeof(XPRMdsofct), tabfct,
     sizeof(tabtyp)/sizeof(XPRMdsotyp), tabtyp,
     sizeof(tabserv)/sizeof(XPRMdsoserv), tabserv
    }; 

In the inverse direction, the following two structures are employed for storing the API information received from Mosel and the inter-module communication interface (IMCI) of mmnl:

static XPRMnifct mm;       /* For storing the Mosel NI function table */
static mmnl_imci mmnl;     /* mmnl API */

Initialization function

The module initialization function performs the initialization of the solver library, communicates the module functionality (interface structures) to Mosel and retrieves the API definitions of the Mosel NI.

DSO_INIT myxprs_init(XPRMnifct nifct, int *interver,int *libver,
  XPRMdsointer **interf)
{
 int r;

 *interver=XPRM_NIVERS;         /* The interface version we are using */
 *libver=XPRM_MKVER(0,0,1);     /* The version of the module: 0.0.1 */
 *interf=&dsointer;             /* Our module interface structure */
 r=XPRSinit(NULL);              /* Initialize the solver */
 if((r!=0)&&(r!=32))
 {
  nifct->dispmsg(NULL,"myxprs: I cannot initialize Xpress Optimizer.\n");
  return 1;
 }
 mm=nifct;                      /* Retrieve the Mosel NI function table */
 return 0;
} 

Implementation of subroutines

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,"myqxprs: 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,"myqxprs: 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 module functions implementing the minimize and maximize subroutines for linear and nonlinear constraints all map to the same function slv_optim.

/**** Linear objective ****/
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,NULL);
}

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,NULL);
}

/**** Nonlinear objective ****/
static int slv_lc_nmaxim(XPRMcontext ctx, void *libctx)
{
 XPRMnlctr obj;

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

static int slv_lc_nminim(XPRMcontext ctx, void *libctx)
{
 XPRMnlctr obj;

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

The function slv_optim first clears any existing solution information (call to slv_clearsol), it then generates the matrix representation of the problem and loads it into the solver (slv_loadmat) and finally starts the actual solving process. After termination of the solver run it retrieves problem status and solution information in order to update the mmnl problem status information via a call to setsolstat.

/**** Optimize a problem: load it, then call the solver ****/
static int slv_optim(XPRMcontext ctx, s_slvctx *slctx, int objsense,
  XPRMlinctr obj, XPRMnlctr nlobj)
{
 int c,i;
 XPRMnlpdata nlpdata;
 s_slvpb *slpb;
 int result;
 double objval;

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

 /* Call function 'loadmat' to generate and load the matrix */
 nlpdata=slv_loadmat(ctx,slctx,obj,nlobj);

 if(nlpdata==NULL)
 {
  mm->dispmsg(ctx,"myqxprs: loadprob failed.\n");
  slpb->saved_ctx=NULL;
  return XPRM_RT_ERROR;
 }

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

 mmnl->setsolstat(ctx,nlpdata,XPRM_PBSOL,0);  /* Solution available for callbacks */
 if((nlpdata->ngents==0) && (nlpdata->nsos==0))
 {
  /* Solve an LP problem */
  c=XPRSlpoptimize(slpb->xpb,"");
  if(c!=0)
  {
   mm->dispmsg(ctx,"myqxprs: optimisation failed.\n");
   slpb->saved_ctx=NULL;
   return XPRM_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,"myqxprs: optimisation failed.\n");
   slpb->saved_ctx=NULL;
   return XPRM_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;
   default:
      result=XPRM_PBOTH;
  }
  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 */
 mmnl->setsolstat(ctx,nlpdata,result,objval);
 slpb->saved_ctx=NULL;
 return 0;
}

The module routine slv_clearsol frees up solution information held in the solver problem interface structures.

/**** Delete solution information ****/
static void slv_clearsol(XPRMcontext ctx, void *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 matrix generation routine invoked by slv_optim looks as follows (given that this implementation works with a quadratic solver any general nonlinear constraints or objectives are rejected). It first calls the mmnl subroutine loadprob to generate the matrix from the problem representation in Mosel. The maxtrix information is then copied into the structures that are expected by the solver's problem input routine (in our case: XPRSloadqcqpglobal). The enumeration employs the mmnl functions buildqtls, freeqtls, and strip.

/**** Load the matrix into the solver ****/
static XPRMnlpdata slv_loadmat(XPRMcontext ctx, s_slvctx *slctx, XPRMlinctr obj,
  XPRMnlctr nlobj)
{
 XPRMnlpdata nlpdata;
 s_slvpb *slpb;
 static char ctcvr[7]={'N','G','L','E','R','1','2'};
 int c,r;
 char pbname[80];
 double *dqe;
 int *mqc1,*mqc2,nb_qp;
 int *qcrows,*qcnqtr,*qcmqc1,*qcmqc2;
 double *qcdqe;
 double objval;

 slpb=SLVCTX2PB(slctx);
 /* Call 'mmnl' function 'loadprob' to generate the matrix */
 nlpdata=NULL;
 c=mmnl->loadprob(ctx,*(slctx->nlctx),&nlpdata,
    NLOPT_COLWJAC|NLOPT_QUADLIN|
    ((slctx->options&OPT_VERBOSE)?NLOPT_VERBOSE:0)|
    ((slctx->options&OPT_LOADNAMES)?NLOPT_NAMES:0),
    obj,nlobj,NULL,NULL,NULL);
 if(c!=0)
  return NULL;

 if(nlpdata->nbqctr<nlpdata->nbnlctr)
 {
  mm->dispmsg(ctx,"myqxprs: Problem contains nonlinear constraints.\n");
  mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
  return NULL;
 }
 if(nlpdata->objtype==1)        /* Objective is nonlinear */
 {
  mm->dispmsg(ctx,"myqxprs: Objective is nonlinear.\n");
  mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
  return NULL;
 }

 /* Define the callback functions for exposing the solution values */
 /* The solution arrays will be stored into our problem object (slpb) */
 nlpdata->fctctx=slpb;
 nlpdata->getsol_v=slv_getsol_v;
 nlpdata->getsol_c=slv_getsol_c;

 /* Xpress Optimizer uses a character code for each variable type */
 for(c=0;c<nlpdata->nbctr;c++)
  nlpdata->ctrtype[c]=ctcvr[(int)(nlpdata->ctrtype[c])];

 if(nlpdata->objtype==2)        /* Quadratic objective */
 {
  XPRMnlqterm qtls;

  r=mmnl->buildqtls(ctx,*(slctx->nlctx),nlpdata->nlobj,1,&nb_qp,&qtls,NULL);
  if((r<0)||(r>=4))
  {
   mmnl->freeqtls(ctx,*(slctx->nlctx),qtls);
   mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
   return NULL;
  }

  if(nb_qp>0)
  {
   if((dqe=malloc((sizeof(double)+sizeof(int)*2)*nb_qp))==NULL)
   {
    mmnl->freeqtls(ctx,*(slctx->nlctx),qtls);
    mm->dispmsg(ctx,"myqxprs: Out of memory error.\n");
    mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
    return NULL;
   }
   mqc1=(int *)(dqe+nb_qp);
   mqc2=mqc1+nb_qp;
   for(c=0;c<nb_qp;c++)
   {
    mqc1[c]=mm->getvarnum(ctx,qtls[c].v1);
    mqc2[c]=mm->getvarnum(ctx,qtls[c].v2);
    dqe[c]=(qtls[c].v1==qtls[c].v2)?qtls[c].coeff*2:qtls[c].coeff;
   }
  }
  else
  {
   mqc1=mqc2=NULL;
   dqe=NULL;
  }
  mmnl->freeqtls(ctx,*(slctx->nlctx),qtls);
 }
 else
 {
  nb_qp=0;
  mqc1=mqc2=NULL;
  dqe=NULL;
 }

 qcrows=qcnqtr=qcmqc1=qcmqc2=NULL;
 qcdqe=NULL;
 if(nlpdata->nbqctr>0)          /* Quadratic constraints */
 {
  XPRMnlqterm qtls;
  int maxterms,curt;
  double rhs;

  qcrows=(int *)malloc(sizeof(int)*2*nlpdata->nbqctr);
  if(qcrows==NULL)
  {
   mm->dispmsg(ctx,"myqxprs: Out of memory error.\n");
   free(dqe);
   mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
   return NULL;
  }
  qcnqtr=qcrows+nlpdata->nbqctr;

  maxterms=0;
  curt=0;
  for(c=0;c<nlpdata->nbqctr;c++)
  {
   qcrows[c]=c;
   r=mmnl->buildqtls(ctx,*(slctx->nlctx),nlpdata->nlctr[c],1,
                                                     &(qcnqtr[c]),&qtls,&rhs);
   if((r<0)||(r>=4))            /* Should not be possible... */
   {
    mmnl->freeqtls(ctx,*(slctx->nlctx),qtls);
    free(qcdqe);
    free(qcmqc2);
    free(qcmqc1);
    free(qcrows);
    free(dqe);
    mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
    return NULL;
   }

   if(qcnqtr[c]+curt>maxterms)
   {
    int *newp1,*newp2;
    double *newp3;
    int newsize;

    newsize=maxterms+ALIGN128(qcnqtr[c]);
    newp1=newp2=NULL;
    newp3=NULL;
    if(((newp1=(int *)realloc(qcmqc1,newsize*sizeof(int)))==NULL)||
       ((newp2=(int *)realloc(qcmqc2,newsize*sizeof(int)))==NULL)||
       ((newp3=(double *)realloc(qcdqe,newsize*sizeof(double)))==NULL))
    {
     mm->dispmsg(ctx,"myqxprs: Out of memory error.\n");
     mmnl->freeqtls(ctx,*(slctx->nlctx),qtls);
     free((newp3!=NULL)?newp3:qcdqe);
     free((newp2!=NULL)?newp2:qcmqc2);
     free((newp1!=NULL)?newp1:qcmqc1);
     free(qcrows);
     free(dqe);
     mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
     return NULL;
    }
    else
    {
     qcmqc1=newp1;
     qcmqc2=newp2;
     qcdqe=newp3;
     maxterms=newsize;
    }
   }

   for(r=0;r<qcnqtr[c];r++,curt++)
   {
    qcmqc1[curt]=mm->getvarnum(ctx,qtls[r].v1);
    qcmqc2[curt]=mm->getvarnum(ctx,qtls[r].v2);
    qcdqe[curt]=(qtls[r].v1==qtls[r].v2)?qtls[r].coeff:qtls[r].coeff/2;
   }
   nlpdata->rhs[c]=-rhs;
   mmnl->freeqtls(ctx,*(slctx->nlctx),qtls);
  }
 }

 /* Load the matrix into the optimizer */
 sprintf(pbname,"xpb%p",slpb);
 r=XPRSloadqcqpglobal(slpb->xpb,pbname,nlpdata->nbvar,
         nlpdata->nbctr,
         nlpdata->ctrtype,nlpdata->rhs,nlpdata->range,nlpdata->grdobj,
         nlpdata->colstart,NULL,nlpdata->cref,
         nlpdata->jac,nlpdata->dlb,nlpdata->dub,
         nb_qp,mqc1,mqc2,dqe,
         nlpdata->nbqctr,qcrows,qcnqtr,qcmqc1,qcmqc2,qcdqe,
         nlpdata->ngents,nlpdata->nsos,nlpdata->qgtype,
         nlpdata->mgcols,nlpdata->mplim,
         nlpdata->qstype,
         nlpdata->msstart,nlpdata->mscols,nlpdata->dref);
 free(qcdqe);
 free(qcmqc2);
 free(qcmqc1);
 free(qcrows);
 free(dqe);

 if(!r)
 {
  /* Set objective constant term */
  c=-1;
  objval=-nlpdata->fixobj;
  XPRSchgobj(slpb->xpb,1,&c,&objval);

 /* Release matrix information - it is no longer required */
 mmnl->strip(ctx,nlpdata,NLSTRIP_ALL);
 return nlpdata;
}

Implementation of mmnl solver interface functions

The following functions implement the mmnl solver interface functions 'getsol_v' and 'getsol_c' that are communicated to mmnl via the NLP problem definition interface structure of type XPRMnlpdata (for further detail see the documentation of loadprob). 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 *mctx, int what, int col)
{
 s_slvpb *slpb;

 slpb=mctx;
 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,"myqxprs: 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,"myqxprs: 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 *mctx, int what, int row)
{
 s_slvpb *slpb;

 slpb=mctx;
 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,"myqxprs: 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,"myqxprs: 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 myqxprs interface for a run ****/
static void *slv_reset(XPRMcontext ctx, void *libctx)
{
 XPRMdsolib dso;
 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,"myqxprs: Out of memory error.\n");
   return NULL;
  }
  memset(slctx,0,sizeof(s_slvctx));

  /* Load interface functions and context of 'mmnl' */
  if(((dso=mm->finddso("mmnl"))==NULL)||
     ((slctx->nlctx=mm->getdsoctx(ctx,dso,(void **)(&(mmnl))))==NULL))
  {
   mm->dispmsg(ctx,"myqxprs: mmnl not found.\n");
   free(slctx);
   return NULL;
  }

  /* Record the problem ID of our problem type */
  mm->gettypeprop(ctx,mm->findtypecode(ctx,"mpproblem.mqxp"),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 those of other example modules (see Mosel NI User Guide), they are listed here for completeness' sake.

/**** 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. 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 QCQP solver example with Xpress Optimizer we have implemented the type handling routines to create, delete, and copy optimization problems (for a full list of type handling routines the reader is referred to the Section 'Table of types' of the Mosel NI Reference Manual). 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,"myqxprs: 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,"myqxprs: I cannot create the problem.\n");
  free(slpb);
  return NULL;
 }
 slpb->slctx=slctx;
 /* Redirect solver messages to the Mosel streams */
 XPRSaddcbmessage(slpb->xpb,slv_cb_output,slpb,0);
 XPRSsetintcontrol(slpb->xpb,XPRS_OUTPUTLOG,1);

 /* Define log callbacks to report program interruption */
 XPRSaddcblplog(slpb->xpb,(void*)slv_cb_stopxprs,slpb,0);
 XPRSaddcbcutlog(slpb->xpb,(void*)slv_cb_stopxprs,slpb,0);
 XPRSaddcbgloballog(slpb->xpb,(void*)slv_cb_stopxprs,slpb,0);
 XPRSaddcbbarlog(slpb->xpb,(void*)slv_cb_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 and it clears any solution information that might have been stored for this problem.

/**** 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;
}

Generating names for matrix entries and matrix output

While developing a solver interface it may be helpful to have the possibility of writing out a matrix representation of the problem held in the solver.

In our example implementation the matrix gets loaded into the solver through the call to optimization, so writing out the problem matrix from the solver needs to take place after this call:

  setparam("myxp_loadnames", true)
  maximize(MyObj)
  writeprob("mymat.lp", "l")

When writing out a matrix for debugging purposes one might expect to be able to match the rows and columns to the Mosel modeling entities via their respective names. By default, Mosel does not generate any names for decision variables or constraints in order to maintain a low memory footprint. This feature needs to be added explicitly into the implementation of the slv_loadmat routine that we have seen earlier in this chapter. After a call to the NI function genmpnames the function slv_loadnames collects the resulting names are into the corresponding data structures that are expected by the solver library (uploading names for rows, columns, and SOS separately in the case of Xpress Optimizer).

/**** Load the matrix into the optimizer ****/
static XPRMnlpdata slv_loadmat(XPRMcontext ctx, s_slvctx *slctx, XPRMlinctr obj,
  XPRMnlctr nlobj)
{
 /* ... load the problem via mmnl routine 'loadprob' ... */

 /* Generate names for the linear part of the problem */
 if(slctx->options&OPT_LOADNAMES)
  mm->genmpnames(ctx,MM_KEEPOBJ,NULL,0);

 /* ... load the problem matrix into the solver ... */

  /* Load names if requested */
  if(slctx->options&OPT_LOADNAMES)
   slv_loadnames(ctx,slpb,nlpdata);

 /* ... release matrix information that is no longer required... */
}

/**** Load names into the optimizer ****/
static void slv_loadnames(XPRMcontext ctx, s_slvpb *slpb, XPRMnlpdata nlpdata)
{
 char *names,*n;
 size_t totlen,totlen2;
 size_t l;
 int c,r;

 totlen=0;
 for(c=0;c<nlpdata->nbvar;c++)
 {
  l=strlen(mm->getmpname(ctx,MM_MPNAM_COL,c));
  totlen+=l+1;
 }
 totlen2=0;
 for(c=0;c<nlpdata->nb_nlctrs;c++)
 {
  l=strlen(nlpdata->names[c]);
  totlen2+=l+1;
 }
 for(;c<nlpdata->nbctr;c++)
 {
  l=strlen(mm->getmpname(ctx,MM_MPNAM_ROW,c-nlpdata->nb_nlctrs));
  totlen2+=l+1;
 }
 if(totlen<totlen2) totlen=totlen2;
 totlen2=0;
 for(c=0;c<nlpdata->nsos;c++)
 {
  l=strlen(mm->getmpname(ctx,MM_MPNAM_SOS,c));
  totlen2+=l+1;
 }
 if(totlen<totlen2) totlen=totlen2;

 if((names=malloc(totlen))==NULL)
  mm->dispmsg(ctx,"myqxprs: Not enough memory for loading the names.\n");
 else
 {
  n=names;
  for(c=0;c<nlpdata->nbvar;c++)
   n+=strlen(strcpy(n,mm->getmpname(ctx,MM_MPNAM_COL,c)))+1;
  if((r=XPRSaddnames(slpb->xpb,2,names,0,nlpdata->nbvar-1))!=0)
   mm->dispmsg(ctx,"myqxprs: Error when executing `%s'.\n","addnames");
  if(!r && (nlpdata->nbctr>0))
  {
   n=names;
   for(c=0;c<nlpdata->nb_nlctrs;c++)
    n+=strlen(strcpy(n,nlpdata->names[c]))+1;
   for(;c<nlpdata->nbctr;c++)
    n+=strlen(strcpy(n,mm->getmpname(ctx,MM_MPNAM_ROW,c-nlpdata->nb_nlctrs)))+1;
   if((r=XPRSaddnames(slpb->xpb,1,names,0,nlpdata->nbctr-1))!=0)
    mm->dispmsg(ctx,"myqxprs: Error when executing `%s'.\n","addnames");
  }
  if(!r && (nlpdata->nsos>0))
  {
   n=names;
   for(c=0;c<nlpdata->nsos;c++)
    n+=strlen(strcpy(n,mm->getmpname(ctx,MM_MPNAM_SOS,c)))+1;
   if((r=XPRSaddnames(slpb->xpb,3,names,0,nlpdata->nsos-1))!=0)
    mm->dispmsg(ctx,"myqxprs: Error when executing `%s'.\n","addnames");
  }
  free(names);
 }
}

The loading of names is triggered from slv_loadmat subject to the presence of the option flag LOADNAMES that is set via a new module parameter myxp_loadnames which is declared via the following entry in the list of parameters ('myxprsparams'):

 {"myxp_loadnames", XPRM_TYP_BOOL|XPRM_CPAR_READ|XPRM_CPAR_WRITE} 

The writeprob routine shown in the Mosel model extract at the beginning of this section needs to be declared in the list of subroutines structure ('tabfct') by adding a line like the following:

  {"writeprob",2103,XPRM_TYP_NOT,2,"ss",slvlc_writepb}

And the actual implementation in the function slv_lc_writepb consists of a call the the solver's matrix output function, along with some error handling such as a check for write access to the specified location:

static int slv_lc_writepb(XPRMcontext ctx, void *libctx)
{
 s_slvpb *slpb;
 int rts;
 char *dname,*options;
 char ename[MM_MAXPATHLEN];

 slpb=SLVCTX2PB((s_slvctx*)libctx);
 dname=MM_POP_REF(ctx);
 options=MM_POP_REF(ctx);
 if((dname!=NULL)&&            /* Make sure that the file can be created */
    (mm->pathcheck(ctx,dname,ename,MM_MAXPATHLEN,MM_RCHK_WRITE|MM_RCHK_IODRV)==0))
 {
  slpb->saved_ctx=ctx;         /* Save current context for callbacks */
  rts=XPRSwriteprob(slpb->xpb,XNLSconvstrto(XNLS_ENC_FNAME,ename,-1,NULL),options);
  slpb->saved_ctx=NULL;
  if(rts)
  {
   mm->dispmsg(ctx,"myqxprs: Error when executing `writeprob'.\n");
   return RT_IOERR;
  }
  else
   return RT_OK;
 }
 else
 {
  mm->dispmsg(ctx,"myqxprs: Cannot write to '%s'.\n",dname!=NULL?dname:"");
  return RT_IOERR;
 }
}