Solve LP, displaying the initial and optimal tableau

Inputs an MPS matrix file and required optimization sense, and proceeds to solve the problem with minim or maxim. The simplex algorithm is interrupted to get its intial basis, and a tableau is requested with a call to function showtab. Once the solution is found, a second call produces the optimal tableau.

Function showtab retrieves the pivot order of the basic variables, along with other problem information, and then constructs (and displays) the tableau row-by-row using the backwards transformation, btran.

Note that tableau should only be used with matrices whose MPS names are no longer than 8 characters.

/***********************************************************************
   Xpress Optimizer Examples
   =========================

   file tableau.c
   ```````````````
   Read in a user-specified LP and solve it, displaying both the  
   initial and optimal tableau. 

   Inputs an MPS matrix file and required optimisation sense, and  
   proceeds to solve the problem with minim or maxim. The simplex 
   algorithm is interrupted to get its intial basis, and a tableau is 
   requested with a call to function showtab. Once the solution is 
   found, a second call produces the optimal tableau. 
   Function showtab retrieves the pivot order of the basic variables, 
   along with other problem information, and then constructs (and  
   displays) the tableau row-by-row using the backwards transformation,  
   btran.

   Arguments:   int argc          number of arguments (2 or 3 expected) 
                char argv[0]      program name 
                char argv[1]      problem name     
                char argv[2]      objective sense ("min" or "max")  
                                    if this argument is not provided the problem will 
                                    be minimised, by default 

   (c) 2017 Fair Isaac Corporation       
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "xprs.h"

XPRSprob probg;

#define MIN 1                        /* Sense values */
#define MAX 0

void showtab(const char *sState,const char *sProblem);
void XPRS_CC optimizermsg(XPRSprob prob, void* data, const char *sMsg,int nLen,int nMsglvl);
void errormsg(const char *sSubName,int nLineno,int nErrorCode);

int main(int argc,char *argv[])
{
   int    nReturn;                   /* Return value of Xpress subroutine */
   int    nOptimizerVersion;         /* Optimizer version number */
   char   *sProblem;                 /* Problem name */
   char   sLogFile[]="tableau.log";  /* Optimizer log file */
   int    nSense = MIN;              /* Optimisation sense: 1 if min, otherwise max */
   int    nIteration;                /* Number of simplex iterations */
   double dObjValue;                 /* Objective function value */
   int nExpiry;

   char banner[256];

   /* Check the number of arguments is correct */
   if (argc < 2 || argc > 3) {
      printf("Usage: tableau <matrix> [min|max]");
      return 1;
   }

   /* Set the problem name */
   sProblem = argv[1];

   /* Set the optimisation sense (default is minim) */
   if (argc == 3) {
      if (strcmp(argv[2],"max") == 0)
         nSense=MAX;
      else if (strcmp(argv[2],"min") != 0) {
         printf("Usage: 'min' or 'max' expected as second argument\n\n");
         return 1;
      }
   }

   /* Initialise Optimizer */
   nReturn=XPRSinit(NULL);
   XPRSgetbanner(banner); printf("%s",banner);
   if (nReturn == 8) return(1);

   nReturn=XPRScreateprob(&probg);
   if (nReturn != 0 && nReturn != 32) errormsg("XPRScreateprob",__LINE__,nReturn);

   /* Delete and define the Optimizer log file */
   remove(sLogFile);
   if (nReturn=XPRSsetlogfile(probg,sLogFile)) errormsg("XPRSsetlogfile",__LINE__,nReturn);

   /* Get and display the Optimizer version number */
   if (nReturn=XPRSgetintcontrol(probg,XPRS_VERSION,&nOptimizerVersion))
   		errormsg("XPRSgetintcontrol",__LINE__,nReturn);

   printf("Xpress Optimiser Subroutine Library Release %.2f\n\n",
      (float)nOptimizerVersion/100);

   /* Tell Optimizer to call optimizermsg whenever a message is output */
   nReturn=XPRSsetcbmessage(probg,optimizermsg,NULL);

   /* Input the matrix */
   if (nReturn=XPRSreadprob(probg,sProblem,"")) errormsg("XPRSreadprob",__LINE__,nReturn);

   /* Turn presolve off as we want to display the initial tableau */
   if (nReturn=XPRSsetintcontrol(probg,XPRS_PRESOLVE,0)) errormsg("XPRSsetintcontrol",__LINE__,nReturn);

   /* Ask Optimizer not to save solutions to the solution file and to read solutions
      from memory instead */
   if (nReturn=XPRSsetintcontrol(probg,XPRS_SOLUTIONFILE,0)) errormsg("XPRSsetintcontrol",__LINE__,nReturn);

   /* Set the number of simplex iterations to zero */
   if (nReturn=XPRSsetintcontrol(probg,XPRS_LPITERLIMIT,0)) errormsg("XPRSsetintcontrol",__LINE__,nReturn);

   /* Perform the first step in the simplex algorithm */
   if (nReturn=(nSense) ? XPRSminim(probg,"") : XPRSmaxim(probg,""))
      errormsg("minim/maxim",__LINE__,nReturn);

   /* Display the initial tableau */
   showtab("Initial",sProblem);

   /* Continue with the simplex algorithm */
   if (nReturn=XPRSsetintcontrol(probg,XPRS_LPITERLIMIT,1000000))
   		errormsg("XPRSsetintcontrol",__LINE__,nReturn);

   if (nReturn=(nSense) ? XPRSminim(probg,"") : XPRSmaxim(probg,""))
      errormsg("minim/maxim",__LINE__,nReturn);

   /* Get and display the objective function value and the iteration count */
   if (nReturn=XPRSgetdblattrib(probg,XPRS_LPOBJVAL,&dObjValue))
   	 	errormsg("XPRSgetdblattrib",__LINE__,nReturn);

   if (nReturn=XPRSgetintattrib(probg,XPRS_SIMPLEXITER,&nIteration))
   		errormsg("XPRSgetintattrib",__LINE__,nReturn);

   printf("M%simal solution with value %g found in %d iterations.\n",
                                            (nSense) ? "in" : "ax",dObjValue,nIteration);

   /* Display the optimal tableau */
   showtab("Optimal",sProblem);

   /* Release Optimizer */
   if (nReturn=XPRSdestroyprob(probg)) errormsg("XPRSdestroyprob",__LINE__,nReturn);
   if (nReturn=XPRSfree()) errormsg("XPRSfree",__LINE__,nReturn);

   return 0;
}


/************************************************************************************\
* Name:         showtab                                                              *
* Purpose:      Display tableau on screen                                            *
* Arguments:    const char *sState       Problem state                               *
*               const char *sProblem     Problem name                                *
* Return Value: None.                                                                *
\************************************************************************************/

void showtab(const char *sState,const char *sProblem)
{
   int    nReturn;               /* Return value of Optimizer routine */
   int    nRow;                  /* Number of rows */
   int    nCol;                  /* Number of columns */
   int    nSprR;
   int    nVector;               /* Number of vectors: rows + cols + 1 (for RHS) */
   int    namlen;                /* Maximum length of row or column name */
   char   *sVectorName;          /* Row or column name */
   int    nElem;                 /* Number of non-zero matrix elements */
   int    *pRowInd;              /* Row indices of the matrix elements */
   int    iRow;                  /* Row index holder */
   double *pMatElem;             /* Values of the matrix elements */
   int    pColStart[2];          /* Start positions of the columns in pMatElem */
   int    *pPivot;               /* Pivot order of the basic variables */
   int    i,j;                   /* Loop counters */
   double *y,*z,d;               /* Variables to form tableau */
   double *x;                    /* Primal solution values */
   double *slack;                /* Slack values */


   /*** Get problem information ***/

   /* Get the number of rows */
   if (nReturn=XPRSgetintattrib(probg,XPRS_ROWS,&nRow)) errormsg("XPRSgetintattrib",__LINE__,nReturn);

   /* Get the number of columns */
   if (nReturn=XPRSgetintattrib(probg,XPRS_COLS,&nCol)) errormsg("XPRSgetintattrib",__LINE__,nReturn);

   /* Get the number of columns */
   if (nReturn=XPRSgetintattrib(probg,XPRS_SPAREROWS,&nSprR)) errormsg("XPRSgetintattrib",__LINE__,nReturn);

   nVector = nRow + nCol + 1;

   /* Allocate memory to the arrays */
   y           = malloc(nRow    * sizeof(double));
   z           = malloc(nVector * sizeof(double));
   pRowInd     = malloc(nRow    * sizeof(int));
   pMatElem    = malloc(nRow    * sizeof(double));
   pPivot      = malloc(nRow    * sizeof(int));
   x           = malloc(nCol    * sizeof(double));
   slack       = malloc(nRow    * sizeof(double));

   /* Get the maximum name length */
   if (nReturn=XPRSgetintattrib(probg, XPRS_NAMELENGTH, &namlen)) errormsg("XPRSgetintattrib", __LINE__, nReturn);
   namlen *= 8;  /* Because name length is returned in 8-byte units */
   /* Allocate somewhere to store individual names */
   sVectorName = malloc((namlen+1)*sizeof(char));

   /* Check that the memory has been successfully allocated */
   if (!y || !z || !pRowInd || !pMatElem || !sVectorName || !pPivot)
     errormsg("malloc",__LINE__,-1);

   /* Retrieve the pivot order of the basic variables */
   if (nReturn=XPRSgetpivotorder(probg,pPivot)) errormsg("XPRSgetpivotorder",__LINE__,nReturn);

   /*** Construct and display the tableau ***/

   /* Printer header of matrix names */
   printf("\n%s tableau of Problem %s.\n",sState,sProblem);
   printf("Note: slacks on G type rows range from -infinity to 0\n\n     ");
   for (i=0; i<nRow; i++) {
      /* Get and print the individual row names */
      if (nReturn=XPRSgetnamelist(probg,1,sVectorName,namlen+1,NULL,i,i)) errormsg("XPRSgetnamelist", __LINE__, nReturn);
      printf(" %.3s   ",sVectorName);
   }
   for (i=0; i<nCol; i++) {
      /* Get and print the individual column names */
      if (nReturn=XPRSgetnamelist(probg,2,sVectorName,namlen+1,NULL,i,i)) errormsg("XPRSgetnamelist", __LINE__, nReturn);
      printf(" %.3s   ",sVectorName);
   }
   printf(" RHS\n");

   /* Get the tableau row-by-row using the backwards transformation, btran */

   /* For each row iRow in turn, calculate z = e_irow * B^-1 * A */
   for (iRow=0; iRow<nRow; iRow++) {

      /* Set up the unit vector e_irow */
      for (i=0; i<nRow; i++) y[i]=0.0;
      y[iRow]=1.0;

      /* y = e_irow * B^-1 */
      if (nReturn=XPRSbtran(probg,y)) errormsg("XPRSbtran",__LINE__,nReturn);

      /* Form z = y * A for each slack column of A */
      for (i=0; i<nRow; i++) z[i]=y[i];

      /* Form z = y * A, for each structural column of A */
      for (j=0; j<nCol; j++) {
         if (nReturn=XPRSgetcols(probg,pColStart,pRowInd,pMatElem,nRow,&nElem,j,j))
            errormsg("XPRSgetcols",__LINE__,nReturn);
         for (d=0.0, i=0; i<nElem; i++)
            d += y[pRowInd[i]] * pMatElem[i];
         z[nRow + j]=d;
      }

      /* Form z for RHS */
      if (nReturn=XPRSgetsol(probg,x,slack,NULL,NULL)) errormsg("XPRSgetsol",__LINE__,nReturn);

      if (pPivot[iRow]>=nRow)
         z[nVector-1]=x[pPivot[iRow]-nRow];
      else
         z[nVector-1]=slack[pPivot[iRow]];

      /* Display single tableau row */
      if (pPivot[iRow]<nRow) { 
        /* Pivot is a row */
        if (nReturn=XPRSgetnamelist(probg,1,sVectorName,namlen+1,NULL,pPivot[iRow],pPivot[iRow])) errormsg("XPRSgetnamelist", __LINE__, nReturn);
      }
      else {
        /* Pivot is a column */
        if (nReturn=XPRSgetnamelist(probg,2,sVectorName,namlen+1,NULL,pPivot[iRow]-nRow-nSprR,pPivot[iRow]-nRow-nSprR)) errormsg("XPRSgetnamelist", __LINE__, nReturn);
      }
      printf("%.3s:",sVectorName);
      for (j=0; j<nVector; j++) {
         if (z[j]>=0.1 || z[j]<=-0.1)
            printf("%5.1f  ",z[j]);
         else
            printf("       ");
      }
      printf("\n");
   }
   printf("\n");

   /* Free memory and close files */
   free(y);
   free(z);
   free(pRowInd);
   free(pMatElem);
   free(pPivot);
   free(sVectorName);
}


/************************************************************************************\
* Name:         optimizermsg                                                              *
* Purpose:      Display Optimizer error messages and warnings.                            *
* Arguments:    const char *sMsg    Message string                                   *
*               int nLen            Message length                                   *
*               int nMsgLvl         Message type                                     *
* Return Value: None.                                                                *
\************************************************************************************/

void XPRS_CC optimizermsg(XPRSprob prob, void* data, const char *sMsg,int nLen,int nMsglvl)
{
   switch (nMsglvl) {

      /* Print Optimizer error messages and warnings */
      case 4:       /* error */
      case 3:       /* warning */
         printf("%*s\n",nLen,sMsg);
         break;

      /* Ignore other messages */
      case 2:       /* dialogue */
      case 1:       /* information */
         break;

      /* Exit and flush buffers */
      default:
         fflush(NULL);
     break;
   }
}


/************************************************************************************\
* Name:         errormsg                                                             *
* Purpose:      Display error information about failed subroutines.                  *
* Arguments:    const char *sSubName    Subroutine name                              *
*               int nLineNo             Line number                                  *
*               int nErrCode            Error code                                   *
* Return Value: None                                                                 *
\************************************************************************************/
void errormsg(const char *sSubName,int nLineNo,int nErrCode)
{
   int nErrNo;      /* Optimizer error number */

   /* Print error message */
   printf("The subroutine %s has failed on line %d\n",sSubName,nLineNo);

   /* Append the error code if it exists */
   if (nErrCode!=-1) printf("with error code %d.\n\n",nErrCode);

   /* Append Optimizer error number, if available */
   if (nErrCode==32) {
      XPRSgetintattrib(probg,XPRS_ERRORCODE,&nErrNo);
      printf("The Optimizer eror number is: %d\n\n",nErrNo);
   }

   /* Free memory close files and exit */
   XPRSdestroyprob(probg);
   XPRSfree();
   exit(nErrCode);
}