mp_qpbandb.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <omp.h>
#include "coiheader.h"
#include "comdecl.h"
 
#ifdef _MSC_VER
#pragma warning (disable : 3280)
#endif
 
/*  QP Model: Min (x-target)*Q*(x-target)/2
    s.t.
              sum( x ) = 1;
              x(i) =l= bin(i);
              sum( bin ) <= maxbin
    where x is a vector of non-negative variables of length NN
          bin is a vector of binary variables of length NN
          target is a vector of +10, and
          Q is a matrix with +1 on the diagonal and 0.1 on the
            first upper and lower bi-diagonal
          maxbin is a constant upper bound on the number of x-components
                that can be nonzero                                   */
 
/* Define Size of the QP model                             */
#define NN 32
#define NVAR   (NN*2)        /* Number of All variables, continuous and binary   */
#define NBIN   NN            /* Number of binary variables                       */
#define MAXBIN 15            /* Max number of nonzero X-variables                */
#define NOBS   60            /* Number of observations used to estimate Q      */
#define NQ     (NN*(NN+1)/2) /* Number of nonzeros in triangle of Q        */
double Q[NQ];
double Mean[NN];
double Target;
int binmap[NN];
int seed = 12359;
 
coiHandle_t *CntVect1;     /* Pointers to CONOPT control vectors, one for each thread.
                             Is global so it can be used in the
                             Solve_Workpool routine                   */
 
void c_log( char *msgt, int code )
/*  This routine has some of the same functionality as the c_log routine in std.c
    but it does not call coiCreate and CoiFree. This is done in the main program. */
{
   FILE *fc;
   int i;
   int c;
 
   if ( code == START )
   {
      i = 0; /* Read the name of the program without extension from stdin */
      while ( (c = getchar()) != EOF && i < MAXLINE ) pname[i++] = (char)c;
      pname[i-1] = '\0'; /* remove eof and new-line chars */
      strcpy(fname,pname); strcat(fname,".lst" );
      fd = fopen(fname,"w");
      strcpy(fname,pname); strcat(fname,".sta" );
      fs = fopen(fname,"w");
   };
   strcpy(fname,pname); strcat(fname,".rc" );
   fc = fopen(fname,"w");
   fprintf(fc,"%s: %s.\n", pname, msgt );
   fclose(fc);
   if ( code != START )
   {
      fclose(fd);
      fclose(fs);
      if ( code == OK )
         exit(0);
      else
         exit(code);
   }
}
 
/* The following global data are general data structures used to manage
   a Branch and Bound process                                         */
 
struct BBNode {
   double *values;/* Solution value or initial values               */
   int    *fixed; /* Binary: 0=fixed zero, 1=fixed 1, -1 not fixed. */
   struct BBNode *prvnode; /* Pointer forward and backward in       */
   struct BBNode *nxtnode; /* double-linked list of active nodes.   */
   double object; /* Objective value (if solved) or lower bound if
                     not yet solved                                 */
   int    node_stat; /* Node Status:
                        0 - created but not solved
                        1 - solved                                  */
   int    solstat;   /* Solution status. Initialized to -1 and
                        redefined when the model has been solved    */
   int    modstat;   /* Model status. Initialized to -1 and
                        redefined when the model has been solved    */
   int    number;    /* Sequense number of the node.                */
   int    depth;     /* Depth of the node in the B&B tree           */
   int    isinteger; /* Is true iff the solution is integer         */
};  /* end struct BBNode                                            */
 
double bintol = 1.0e-7;
int totnodes;  /* Total number of nodes in use                      */
int freenodes; /* Number of active nodes                            */
int actnodes;  /* Number of free nodes in a free list               */
int numsolve;  /* Counter for finished Solves.                      */
/*
   Data for best solution
*/
int intfound;     /* Global flag indicating whether we already have an
                     integer feasible solution.                     */
int numintfound;  /* Number of integer solutions found              */
double bestint;   /* Objective value for the best integer solution
                     found so far.                                  */
double solution[NVAR]; /* The solution corresponding to BestInt          */
/*
   Data for managing the Branch and Bound Tree
*/
struct BBNode *actnode;  /* Pointer to first node in a list of active
                            untested noded                          */
struct BBNode *freenode; /* Pointer to first node in a freelist     */
/*
   Data for managing a WorkPool with the nodes selected for the next
   set of parallel Solves
*/
#define maxmaxpool 64 /* Allocation size for WorkPool                */
int maxpool;         /* Largest size of WorkPool selected by user   */
int workpool;        /* Size of the actual workpool                 */
struct BBNode *work_node[maxmaxpool]; /* Pointers to vector of selected nodes
                              in the WorkPool                       */
int parallel;        /* If set then run the COI_Solve loop over
                        the nodes in the WorkPool in parallel       */
double stime;        /* Accumulated time in Solve Loop              */
 
/* --------------START OF GENERAL BRANCH AND BOUND ROUTINES           */
/* The following routines are general Branch and Bound routines
   based on the global data above                                     */
 
void list_node( struct BBNode *node)
{
 
   if ( node )
      printf("list_node: number= %d node_stat= %d.\n", node->number,node->node_stat);
   else
      printf("list_node: node is NULL.\n");
}
 
void list_nodes()
{
   struct BBNode *node;
   int n = 0;
 
   node = actnode;
   while ( node )
   {
      printf("list_nodes: node number= %d  node_stat= %d  object= %f\n", node->number, node->node_stat, node->object);
      node = node->nxtnode;
      n++;
   }
   printf("list_nodes: Total= %d.\n", n);
}
 
void remove_activenode(struct BBNode *node)
{
/* Remove a node from the list of active nodes.                       */
 
   if ( node->prvnode )
      node->prvnode->nxtnode = node->nxtnode;
   else
      actnode = node->nxtnode;
   if ( node->nxtnode)
      node->nxtnode->prvnode = node->prvnode;
   actnodes--;
}
 
void return_freenode(struct BBNode *node)
{
/* Link the node into the free list. The free list does not use
   prvnode                                                            */
 
   node->nxtnode = freenode;
   freenode      = node;
   freenodes++;
}
 
struct BBNode *get_freenode()
{
/* Get a free node to be used in the list of active nodes and
   initialize the status.                                             */
 
   struct BBNode *node;
   if ( freenode )
   {
      node      = freenode;
      freenode  = freenode->nxtnode;
      freenodes--;
   }
   else
   {
      node = (struct BBNode *)malloc( sizeof(struct BBNode) );
      node->nxtnode = NULL;
      node->prvnode = NULL;
      node->values  = (double *)malloc( NVAR*sizeof(double) );
      node->fixed   = (int    *)malloc( NBIN*sizeof(int) );
      totnodes++;
      node->number = totnodes;
   }
   node->node_stat = 0;
   node->solstat   =-1;
   node->modstat   =-1;
   return node;
}
 
void add_activenode (struct BBNode *node)
{
/* Add a node to the front of the list of Active nodes                */
 
   actnodes++;
 
   if (actnode)
      actnode->prvnode = node;
   node->nxtnode = actnode;
   node->prvnode = NULL;
   actnode = node;
}
 
void add2workpool( struct BBNode *node )
{
/*  Add a node to the Workpool list                                   */
 
   work_node[workpool] = node;
   workpool++;
}
 
int select_var(struct BBNode *node)
{
/* Selects variable to branch on, counted in the space of binary
   variables. Currently there is only one rule, namely to select based
   on max infeasibility. Other Variable selection rules could be
   implemented in this routine.                                       */
 
   int ib = -1, i,j;
   double dist, maxdist = 0;
 
   for ( i=0; i<NBIN; i++ )
   {
      j = binmap[i];
      if ( node->values[j] < 1.0-node->values[j] )
         dist = node->values[j];
      else
         dist = 1.0-node->values[j];
      if ( dist > maxdist )
      {
         maxdist = dist;
         ib      = i;
      }
   }
   return ib;
}
 
void create_workpool()
{
/* Create a workpool of up to MaxPool nodes to be processed in
   parallel.
   The nodes are removed from the active list.                        */
 
   struct BBNode *node, *curr, *next;
   double cutoff;
   int ib;    /* Branching variable                                   */
   int i;
 
/*   printf("Starting create_workpool. maxpool=%d\n", maxpool);
   list_nodes(); */
   node = NULL;
   workpool = 0;
   while ( workpool < maxpool )
   {
      cutoff = 1.e30;
      curr   = actnode;
/*
   Run through the linked list of active nodes and find the best.
   We use the objective but other criteria could be used
*/
      node = NULL;
      while ( curr )
      {
         next = curr->nxtnode;
         if ( curr->object < cutoff )
         {
            node   = curr;
            cutoff =  curr->object;
         }
         curr = next;
      }
      if ( node == NULL ) break;
      printf("Selected node %d with object %f and node_stat %d\n",node->number, node->object, node->node_stat);
      if ( node->node_stat == 0 )
      {
/*
   The node has been created but not solved.
   Take it immediately: remove it from the active list and add it
   to the work_node vector
*/
         remove_activenode(node);
         add2workpool(node);
      }
      else
      {
/*
   The node has been solved. Split it in two and take one or both
   depending on the space left in the workpool
 
   Get a new node and, select a branching variable, and copy the
   old node into the new node.
*/
         next = get_freenode();
         ib = select_var( node );
         for ( i=0; i<NVAR; i++ )
            next->values[i] = node->values[i];
         for ( i=0; i<NBIN; i++ )
            next->fixed[i]  = node->fixed[i];
         next->object = node->object;
/*
   Make the two node unsolved and fix the branching variable down and up
   respectively.
*/
         node->node_stat = 0;
         node->solstat   = -1;
         node->modstat   = -1;
         node->depth++;
         node->fixed[ib] = 0; /*  branch down */
 
         next->node_stat = 0;
         next->solstat   = -1;
         next->modstat   = -1;
         next->depth     = node->depth;
         next->fixed[ib] = 1;   /* Branch up  */
/*
   Add the first node to the workpool
*/
         remove_activenode(node);
         add2workpool(node);
/*
   Add the second node to the workpool or to the pool of active
   nodes, depending on space in the workpool
*/
         if ( workpool < maxpool )
            add2workpool(next);
         else
            add_activenode(next);
      }
   }
}
 
void solve_workpool()
{
/* Solve the models in the WorkPool and integrate the solutions into
   the active list of nodes. Update the best solution if one is found.*/
 
   int nw;      /* index for active node in the WorkPool              */
   int i;
   int thread;
   int COI_Error;       /* Local Error counter  */
   coiHandle_t CntVect; /* Local control vector */
   struct BBNode *node, *curr, *next;
   int newint;
   double time0;
 
   time0 = omp_get_wtime();
   if ( workpool <= 1 || 0 == parallel )
   {
/*  This is a serial version                         */
      thread = 0;
      for (nw = 0; nw < workpool; ++nw )
      {
         node = work_node[nw];
/*         printf("Starting to solve node %d located at %p.\n",node->number, node); */
         CntVect = CntVect1[thread];
         COI_Error = COIDEF_UsrMem( CntVect, (void *)node ); /* Pass the current node through Usrmem */
         if ( COI_Error )
            c_log("**** Fatal Error while loading CONOPT Callback routines.", -1 );
         COI_Error = COI_Solve( CntVect );
         if ( COI_Error )
            c_log("**** Fatal Error while solving.", -1);
      }
   }
   else
/* And this is the corresponding Parallel version of the SOLVE loop.  */
   {
      COI_Error = 0;
#pragma omp parallel for default(shared) private(CntVect,thread,node) schedule(dynamic) reduction(+:COI_Error)
      for (nw = 0; nw < workpool; ++nw )
      {
         thread = omp_get_thread_num();
         node = work_node[nw];
         printf("Solving node %d with thread %d.\n", node->number, thread);
         CntVect = CntVect1[thread];
         COI_Error = COIDEF_UsrMem( CntVect, (void *)node ); /* Pass the current node through Usrmem */
         if ( !COI_Error )
            COI_Error = COI_Solve( CntVect );
      }
      if ( COI_Error ) /* We cannot exit inside a parallel loop so we wait */
         c_log("**** Fatal Error inside parallel loop.", -1);
   }
/*
   Endif */
   stime = stime + omp_get_wtime() - time0;
   numsolve += workpool;
/*
   Process the solutions. The loop is not parallelized so the order here
   is fixed, independent of the order in which the nodes finished
   solving.
*/
   newint = 0;  /* No new integer solution in this pass               */
   for (nw = 0; nw < workpool; ++nw )
   {
      node = work_node[nw];
      printf("Node= %d Solstat= %d Modstat= %d Object= %f Isinteger = %d\n",
         node->number, node->solstat, node->modstat, node->object, node->isinteger);
      if ( node->solstat != 1 || node->modstat > 2 )
      {
/*  The node was not solved to local optimality. Fathom.              */
         printf("Not optimal -- fathom.\n");
         return_freenode(node);
      }
      else if ( intfound )
      {
         if ( node->object >= bestint )
         {
/*   A better solution has been found. Fathom.                        */
            printf("Worse than optimal -- fathom.\n");
            return_freenode(node);
         }
         else
         {
            if ( node->isinteger )
            {
               printf("A: Integer solution stored.\n");
               numintfound++;
               newint   = 1;
               bestint  = node->object;
               for ( i = 0; i < NVAR; i++ )
                  solution[i] = node->values[i];
               return_freenode(node);
            }
            else
            {
               printf("Non-integer. kept.\n");
               add_activenode(node);
            }
         }
      }
      else
/*
   An Integer solution has not yet been found. We cannot fathom it.
*/
      {
         if ( node->isinteger )
         {
            printf("B: Integer solution stored.\n");
            numintfound++;
            intfound = 1;
            newint   = 1;
            bestint  = node->object;
            for ( i = 0; i < NVAR; i++ )
               solution[i] = node->values[i];
            return_freenode(node);
         }
         else
         {
            printf("Non-integer. kept.\n");
            add_activenode(node);
         }
      }
   }
/*
   If we did find a new integer solution then scan the list of active
   nodes to check if some of them should be fathomed.
*/
   if ( newint )
   {
      curr = actnode;
      while (curr)
      {
         next = curr->nxtnode;
         if ( curr->object > bestint )
         {
            printf("Node %d with object= %f fathomed by new objective.\n", curr->number, curr->object);
            remove_activenode(curr);
            return_freenode(curr);
         }
         curr = next;
      }
   }
}
 
void bandb_solve()
{
   int i;
   struct BBNode *node;
/*
  Initialize the first node for B&B tree
*/
   numsolve    = 0;
   numintfound = 0;
   freenode    = NULL;
   actnode     = NULL;
   intfound    = 0;
   totnodes    = 0;   /* total number of nodes in use                 */
   actnodes    = 0;   /* number of active nodes                       */
   freenodes   = 0;   /* number of free nodes in a free list          */
   stime       = 0.0; /* reset solve loop timer                       */
 
   node = get_freenode();/* Create a root node with no binaries fixed    */
   node->depth = 0;
   for ( i = 0; i < NVAR; i++ )
      node->values[i] = 0.0;
   for ( i = 0; i < NBIN; i++ )
      node->fixed[i] = -1;
   node->object = 1.e29; /* We need a value for the first selection
                            before the node has been solved           */
   add_activenode(node);
   while ( actnodes > 0 )
   {
      printf("Creating WorkPool. Active nodes= %d\n",actnodes);
      create_workpool();
      printf("WorkPool with %d nodes selected. Remaining Active nodes= %d.\n",workpool,actnodes);
      if ( workpool > 0 )
         solve_workpool();
      else if ( workpool == 0 && actnodes > 0 )
         c_log("No nodes found but actnodes is positive", -1);
   }
 
   if ( intfound )
   {
      printf("\nAn integer solution with objective %f has been found.\n\n",bestint);
      printf("Solution values:\n");
      for ( ; i < NVAR; i++ )
         printf("  %d   %f  \n", i, solution[i]);
   }
   else
      printf("No integer solution has been found\n");
   printf("Number of Solves:            %d\n", numsolve);
   printf("Number of Integer Solutions: %d\n\n", numintfound);
 
   while( freenode != NULL )
   {
       struct BBNode* n = freenode->nxtnode;
       free(freenode->values);
       free(freenode->fixed);
       free(freenode);
       freenode = n;
   }
}
/* --------------END OF GENERAL BRANCH AND BOUND ROUTINES             */
float rndx()
{
/* Defines a pseudo random number between 0 and 1                     */
 
   seed = (seed*1027+25)%1048576;
   return ((float)seed)/((float)1048576);
}
 
void defdata()
{
/*  Define values Target, Mean, and Q                                 */
 
   int i, j, j1, k, l;
   double Sum;
   double *Xobs;
 
   Target = 10.0;
   Xobs = (double *)malloc(NN*NOBS*sizeof(double));
 
   for (i=0; i< NN; i++)
      Mean[i] = Target + (rndx()*4.0-2.0);
/*
  Create a positive definite matrix computed as the sum of some
  outer products                                                      */
   j = 0;
   for (i=0; i< NQ; i++)
      Q[i] = 0;
 
   for (i=0; i< NN; i++)
   {
      Sum = 0.0;
      j1 = j;
      for (k=0; k < NOBS; k++)
      {
         Xobs[j] = 100.0*rndx();
         Sum     = Sum + Xobs[j];
         j++;
      }
      Sum = Sum / NOBS;
      j = j1;
      for (k=0; k < NOBS; k++)
      {
         Xobs[j] = Xobs[j] - Sum;
         j++;
      }
   }
/*  Define the Q-matrix as the sum deviations and store it as a lower
    triangular matrix stored column-wise in this order
    0
    1    NN
    2    NN+1
    .
    NN-1 ..                                                           */
 
   l = 0;
   for (i=0; i< NN; i++)
   {
      for (j=i; j < NN; j++)
      {
         Sum = 0.0;
         for (k=0; k < NOBS; k++)
            Sum = Sum + Xobs[i*NN+k] * Xobs[j*NN+k];
         Q[l] = Sum / (NOBS-1);
         l++;
      }
   }
   free(Xobs);
   for ( i=0; i<NN; i++ )
      binmap[i] = NN+i;
 
}

int COI_CALLCONV QP_ReadMatrix( double LOWER[], double CURR[], double UPPER[], int VSTA[], int TYPE[], double RHS[],
                            int ESTA[], int COLSTA[], int ROWNO[], double VALUE[],
                            int NLFLAG[], int NUMVAR, int NUMCON, int NUMNZ, void* USRMEM )
{
  int i, j;
  struct BBNode *node;
  node = (struct BBNode *) USRMEM;
 /*                                                               */
 /*   Information about Variables:                                */
 /*   Default: Lower = -Inf, Curr = 0, and Upper = +inf.          */
 /*   Default: the status information in Vsta is not used.        */
 /*                                                               */
 /*   Lower bound = 0 on all variables:                           */
 /*                                                               */
   for ( i=0; i<NN; i++ )
   {
      LOWER[i] = 0.0 ;
      LOWER[NN+i] = 0.0;
      UPPER[NN+i] = 1.0;
   }
/* Update the current solution and the bounds with information
   for the current node                                           */
   for ( i=0; i<NVAR; i++ )
   {
      CURR[i] = node->values[i];
   }
   for ( i=0; i<NBIN; i++ )
   {
      j = binmap[i];
      if ( node->fixed[i] == 0 )
      {
         LOWER[j] = 0.0;
         UPPER[j] = 0.0;
      }
      else if ( node->fixed[i] == 1 )
      {
         LOWER[j] = 1.0;
         UPPER[j] = 1.0;
      }
   }
 /*                                                               */
 /*   Information about Constraints:                              */
 /*   Default: Rhs = 0                                            */
 /*   Default: the status information in Esta and the function    */
 /*            value in FV are not used.                          */
 /*   Default: Type: There is no default.                         */
 /*            0 = Equality,                                      */
 /*            1 = Greater than or equal,                         */
 /*            2 = Less than or equal,                            */
 /*            3 = Non binding.                                   */
 /*                                                               */
 /*   Constraint 0 (Objective)                                    */
 /*   Rhs = 0.0 and type Non binding                              */
 /*                                                               */
   TYPE[0] = 3    ;
 /*                                                               */
 /*   Constraint 1 (X Sum to 1)                                   */
 /*   Rhs = 1 and type Equality                                   */
 /*                                                               */
   RHS[1]  = 1.0 ;
   TYPE[1] = 0   ;
 /*                                                               */
 /*   Constraint 2 (Sum X*Mean ) = Target                         */
 /*   Rhs = 1 and type Equality                                   */
 /*                                                               */
   RHS[2]  = Target;
   TYPE[2] = 0   ;
 /*                                                               */
 /*   Constraint 3 (Bin Sum to le Maxbin)                         */
 /*   Rhs = 1 and type Less than or Equal                         */
 /*                                                               */
   RHS[3]  = MAXBIN ;
   TYPE[3] = 2   ;
 /*                                                               */
 /*   Constraint group 4 (X le Bin or X - Bin le 0)               */
 /*   Rhs = 0 and type Less than or Equal                         */
 /*                                                               */
   for ( i=0; i<NN; i++ )
   {
      RHS[4+i] = 0.0 ;
      TYPE[4+i] = 2;
   }
 
 /*                                                               */
 /*   Information about the Jacobian. We use the standard method with */
 /*   Rowno, Value, Nlflag and Colsta and we do not use Colno.    */
 /*                                                               */
 /*   Colsta = Start of column indices (No Defaults):             */
 /*   Rowno  = Row indices                                        */
 /*   Value  = Value of derivative (by default only linear        */
 /*                                 derivatives are used)         */
 /*   Nlflag = 0 for linear and 1 for nonlinear derivative        */
 /*            (not needed for completely linear models)          */
 /*                                                               */
   j = 0;  /* counter for current nonzero   */
   for ( i=0; i<NN; i++ ) {  /* X-variables */
      COLSTA[i] = j;
      ROWNO[j]  = 0;  /* Row 0: Objective   */
      NLFLAG[j] = 1;
      j++;
      ROWNO[j]  = 1;  /* Row 1: Sum(X) le 1     */
      VALUE[j]  = 1.0;
      NLFLAG[j] = 0;
      j++;
      ROWNO[j]  = 2;  /* Row 2: Sum(X*Mean) le Target    */
      VALUE[j]  = Mean[i];
      NLFLAG[j] = 0;
      j++;
      ROWNO[j]  = 4+i;  /* Row group 4: X - Bin le 0 */
      VALUE[j]  = 1.0;
      NLFLAG[j] = 0;
      j++;
      }
   for ( i=0; i<NN; i++ ) {  /* Bin-variables */
      COLSTA[NN+i] = j;
      ROWNO[j]  = 3;  /* Row 3: Sum(Bin) le Maxbin     */
      VALUE[j]  = 1.0;
      NLFLAG[j] = 0;
      j++;
      ROWNO[j]  = 4+i;  /* Row group 4: X - Bin le 0 */
      VALUE[j]  = -1.0;
      NLFLAG[j] = 0;
      j++;
      }
   COLSTA[2*NN] = j;
   return 0;
}
 

int COI_CALLCONV QP_FDEval( const double X[], double* G, double JAC[], int ROWNO, const int JCNM[], int MODE,
                          int IGNERR, int* ERRCNT, int NUMVAR, int NUMJAC, int THREAD, void* USRMEM )
{
   int i, j, k;
   double sum;
 /*                                                                   */
 /*   Row 0: the objective function is nonlinear                      */
 /*                                                                   */
 
   if ( ROWNO == 0 )
   {
 /*                                                                   */
 /*   Mode = 1 or 3. Function value: G = x*Q*x                        */
 /*                                                                   */
      if ( MODE == 1 || MODE == 3 )
      {
         k = 0;
         sum = 0.0;
         for (i=0; i< NN; i++)
         {
            sum += X[i]*Q[k]*X[i];
            k++;
            for (j=i+1; j < NN; j++)
            {
               sum += 2*X[i]*Q[k]*X[j];
               k++;
            }
         }
         *G = sum / 2;
      }
 /*                                                                   */
 /*   Mode = 2 or 3: Derivative values:                               */
 /*                                                                   */
      if ( MODE == 2 || MODE == 3 )
      {
         for ( i=0; i<NN; i++ )
            JAC[i] = 0.0;
         k = 0;
         for (i=0; i < NN; i++)
         {
            JAC[i] += Q[k]*X[i];
            k++;
            for (j=i+1; j < NN; j++)
            {
               JAC[i] += Q[k]*X[j];
               JAC[j] += Q[k]*X[i];
               k++;
            }
         }
      }
   }
 /*                                                                   */
 /*   Other rows: These rows are linear and will not be called.       */
 /*                                                                   */
 
   return 0;
}
 

int COI_CALLCONV QP_2DDir( const double X[], const double DX[], double D2G[], int ROWNO,
                const int JACNUM[], int* NODRV, int NUMVAR, int NUMJAC, int THREAD, void* USRMEM )
{
int i, j, k;
/*                                                                    */
/*   Is only called for Irow = 0, the nonlinear objective             */
/*   Return H*Dx = Q*Dx in D2G                                        */
/*                                                                    */
   if ( ROWNO == 0 )
   {
      for ( i=0; i<NN; i++ )
         D2G[i] = 0.0;
      k = 0;
      for (i=0; i< NN; i++)
      {
         D2G[i] += Q[k]*DX[i];
         k++;
         for (j=i+1; j < NN; j++)
         {
            D2G[i] += Q[k]*DX[j];
            D2G[j] += Q[k]*DX[i];
            k++;
         }
      }
   }
   return 0;
}
 

int COI_CALLCONV QP_2DLagrStr( int HSRW[], int HSCL[], int* NODRV,
                               int NUMVAR, int NUMCON, int NHESS, void* USRMEM )
{
   int i, j, k;
 
   k = 0;
   for (i=0; i< NN; i++)
      for (j=i; j < NN; j++)
      {
         HSRW[k] = j;
         HSCL[k] = i;
         k++;
      }
   return 0;
}
 

int COI_CALLCONV QP_2DLagrVal( const double X[], const double U[], const int HSRW[], const int HSCL[],
                double HSVL[], int* NODRV, int NUMVAR, int NUMCON, int NHESS, void* USRMEM )
{
   int k;
 
   for ( k=0; k<NQ; k++ )
      HSVL[k] = Q[k]*U[0];
   return 0;
}
 
int COI_CALLCONV BB_Message( int SMSG, int DMSG, int NMSG, char* MSGV[], void* USRMEM  )
{
/* Simple implementation that does nothing to avoid mixing output from
   multiple threads.                                                  */
   return 0;
}
 
int COI_CALLCONV BB_ErrMsg( int ROWNO, int COLNO, int POSNO, const char* MSG, void* USRMEM )
{
/* Simple implementation that does nothing to avoid mixing output from
   multiple threads.                                                  */
   return 0;
}
 
int COI_CALLCONV BB_Status( int MODSTA, int SOLSTA, int ITER, double OBJVAL, void* USRMEM )
{
/*
   Simple implementation in which we just store the status and
   objective value in the current node.
*/
   struct BBNode *node;
 
   node = (struct BBNode *) USRMEM;
 
   stacalls++;
   node->object  = OBJVAL;
   node->modstat = MODSTA;
   node->solstat = SOLSTA;
   node->node_stat = 1;
   return 0;
}
 
int COI_CALLCONV BB_Solution( const double XVAL[], const double XMAR[], const int XBAS[], const int XSTA[],
                          const double YVAL[], const double YMAR[], const int YBAS[], const int YSTA[],
                          int NUMVAR, int NUMCON, void* USRMEM )
{
/*
   Store the solution in the Values field of the solved node.
   Also check for integrality. If done here it is inside the
   parallel loop.                                                     */
   struct BBNode *node;
   int i;
   int j;
   int numintvar=0;
   double dist;
 
   node = (struct BBNode *) USRMEM;
   for ( i=0; i<NUMVAR; i++ )
      node->values[i] = XVAL[i];
   for ( i=0; i < NBIN; i++ )
   {
      j    = binmap[i];
      if ( node->values[j] < 1.0-node->values[j] )
         dist = node->values[j];
      else
         dist = 1.0-node->values[j];
      if ( dist <= bintol ) numintvar++;
   }
   node->isinteger = (numintvar==NBIN);
   solcalls++;
 
   return 0;
}

int main(int argc, char** argv)
{
   coiHandle_t CntVect; /* the global CntVect is overwritten in this model */
   int maxthread;       /* Maximum number of threads                       */
   char msg[256];
 
   int thread;
   double  time0;
   double  time1;
   double  time8S,  time16S,  time32S;
   double  time8P,  time16P,  time32P;
   double stime1;
   double stime8S, stime16S, stime32S;
   double stime8P, stime16P, stime32P;
   int numsolve1;
   int numsolve8S, numsolve16S, numsolve32S;
   int numsolve8P, numsolve16P, numsolve32P;
 
   c_log( "Starting to execute", START );
   maxthread = omp_get_max_threads();
   printf("Initial maxthread= %d\n",maxthread);
   if ( maxthread < 4 ) {
      maxthread = 4;
      omp_set_num_threads(maxthread);
   }
   printf("Revised maxthread= %d\n",maxthread);
 
/*
   Allocate and initialize the Q matrix
 */
   defdata();
/*
   We cannot use COIDEF_IniC since we need several control vectors.
   Below we initialize with COIDEF_Ini, COIDEF_Base, and COIDEF_C
   for the control record for each block.
*/
   CntVect1 = (coiHandle_t *) malloc( maxthread*sizeof(coiHandle_t) );
 /*
   Tell CONOPT about the sizes in the model
 */
   for ( thread=0; thread<maxthread; thread++ )
   {
      CntVect1[thread] = NULL;
      coiCreate(&CntVect);
      if (NULL==CntVect) {
         printf("Could not create Conopt object: %s\n", msg);
         exit(1);
      }
      CntVect1[thread] = CntVect;
      COI_Error +=  COIDEF_NumVar    ( CntVect, 2*NN );       /* 2*NN variables                    */
      COI_Error +=  COIDEF_NumCon    ( CntVect, 4+NN );       /* 3+NN constraints                  */
      COI_Error +=  COIDEF_NumNz     ( CntVect, 6*NN );       /* 5*NN nonzeros in the Jacobian     */
      COI_Error +=  COIDEF_NumNlNz   ( CntVect, NN   );       /* NN of which are nonlinear         */
      COI_Error +=  COIDEF_NumHess   ( CntVect, NQ   );       /* NQ nonzero hessian elements       */
      COI_Error +=  COIDEF_OptDir    ( CntVect, -1   );       /* Minimize                          */
      COI_Error +=  COIDEF_ObjCon    ( CntVect, 0    );       /* Objective is constraint 0         */
      COI_Error +=  COIDEF_DebugFV   ( CntVect, 0    );       /* 0 means no debugging, 1 each iter */
      COI_Error +=  COIDEF_StdOut    ( CntVect, 0    );       /* 1 means Allow output to StdOut    */
 /*
   Register the necessary callback routines with CONOPT
 */
      COI_Error +=  COIDEF_Message   ( CntVect, &BB_Message  );  /* Register the callback Message     */
      COI_Error +=  COIDEF_ErrMsg    ( CntVect, &BB_ErrMsg   );  /* Register the callback ErrMsg      */
      COI_Error +=  COIDEF_Status    ( CntVect, &BB_Status   );  /* Register the callback Status      */
      COI_Error +=  COIDEF_Solution  ( CntVect, &BB_Solution );  /* Register the callback Solution    */
      COI_Error +=  COIDEF_ReadMatrix( CntVect, &QP_ReadMatrix);  /* Register the callback ReadMatrix  */
      COI_Error +=  COIDEF_FDEval    ( CntVect, &QP_FDEval);      /* Register the callback FDEval      */
      COI_Error +=  COIDEF_2DDir     ( CntVect, &QP_2DDir );      /* Register the callback 2DDir       */
      COI_Error +=  COIDEF_2DLagrStr ( CntVect, &QP_2DLagrStr);   /* Register the callback 2DLagrStr   */
      COI_Error +=  COIDEF_2DLagrVal ( CntVect, &QP_2DLagrVal);   /* Register the callback 2DLagrVal   */
 
#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
      COI_Error +=  COIDEF_License   ( CntVect, LICENSE_INT_1, LICENSE_INT_2, LICENSE_INT_3, LICENSE_TEXT);
#endif
   }
 
   if ( COI_Error ) {
      printf("Skipping COI_Solve due to setup errors. COI_Error = %d\n",COI_Error);
      c_log( "Skipping Solve due to setup errors", COI_Error); }
/*
   Solve the model using sequential processing with 1 node in pool
*/
   parallel = 0;
   maxpool  = 1;
   time0 = omp_get_wtime();
   bandb_solve();
   time1 = omp_get_wtime() - time0;
   stime1 = stime;
   numsolve1 = numsolve;
   fprintf(fd,"Finished sequential with maxpool= %d\n",maxpool);
/*
   Solve the model using sequential processing with 8 nodes in pool
*/
   parallel = 0;
   maxpool  = 8;
   time0 = omp_get_wtime();
   bandb_solve();
   time8S = omp_get_wtime() - time0;
   stime8S = stime;
   numsolve8S = numsolve;
   fprintf(fd,"Finished sequential with maxpool= %d\n",maxpool);
/*
   Solve the model using sequential processing with 16 nodes in pool
*/
   parallel = 0;
   maxpool  = 16;
   time0 = omp_get_wtime();
   bandb_solve();
   time16S = omp_get_wtime() - time0;
   stime16S = stime;
   numsolve16S = numsolve;
   fprintf(fd,"Finished sequential with maxpool= %d\n",maxpool);
/*
   Solve the model using sequential processing with 32 nodes in pool
*/
   parallel = 0;
   maxpool  = 32;
   time0 = omp_get_wtime();
   bandb_solve();
   time32S = omp_get_wtime() - time0;
   stime32S = stime;
   numsolve32S = numsolve;
   fprintf(fd,"Finished sequential with maxpool= %d\n",maxpool);
/*
   Solve the model using sequential processing with 8 nodes in pool
*/
   parallel = 1;
   maxpool  = 8;
   time0 = omp_get_wtime();
   bandb_solve();
   time8P = omp_get_wtime() - time0;
   stime8P = stime;
   numsolve8P = numsolve;
   fprintf(fd,"Finished parallel with maxpool= %d\n",maxpool);
/*
   Solve the model using sequential processing with 16 nodes in pool
*/
   parallel = 1;
   maxpool  = 16;
   time0 = omp_get_wtime();
   bandb_solve();
   time16P = omp_get_wtime() - time0;
   stime16P = stime;
   numsolve16P = numsolve;
   fprintf(fd,"Finished parallel with maxpool= %d\n",maxpool);
/*
   Solve the model using sequential processing with 32 nodes in pool
*/
   parallel = 1;
   maxpool  = 32;
   time0 = omp_get_wtime();
   bandb_solve();
   time32P = omp_get_wtime() - time0;
   stime32P = stime;
   numsolve32P = numsolve;
   fprintf(fd,"Finished parallel with maxpool= %d\n",maxpool);
 
   fprintf(fd,"Number of Binary variables in model= %d\n",NN);
   fprintf(fd,"Number of threads in parallel experiments= %d\n\n",maxthread);
   fprintf(fd,"Loop type  MaxPool  Total time     %%      Solve Time      %%    Number of Solves    Parallel Speedup\n\n");
   fprintf(fd,"Sequential     1    %9.3f    %5.1f    %9.3f     %5.1f          %d\n",
               time1 ,              100.0, stime1 ,                100.0, numsolve1);
   fprintf(fd,"Sequential     8    %9.3f    %5.1f    %9.3f     %5.1f          %d\n",
               time8S, time8S/time1*100.0, stime8S, stime8S/stime1*100.0, numsolve8S);
   fprintf(fd,"Sequential    16    %9.3f    %5.1f    %9.3f     %5.1f          %d\n",
               time16S, time16S/time1*100.0, stime16S, stime16S/stime1*100.0, numsolve16S);
   fprintf(fd,"Sequential    32    %9.3f    %5.1f    %9.3f     %5.1f          %d\n",
               time32S, time32S/time1*100.0, stime32S, stime32S/stime1*100.0, numsolve32S);
   fprintf(fd,"Parallel       8    %9.3f    %5.1f    %9.3f     %5.1f          %d        %10.2f\n",
               time8P, time8P/time1*100.0, stime8P, stime8P/stime1*100.0, numsolve8P, time8S/time8P);
   fprintf(fd,"Parallel      16    %9.3f    %5.1f    %9.3f     %5.1f          %d        %10.2f\n",
               time16P, time16P/time1*100.0, stime16P, stime16P/stime1*100.0, numsolve16P, time16S/time16P);
   fprintf(fd,"Parallel      32    %9.3f    %5.1f    %9.3f     %5.1f          %d        %10.2f\n",
               time32P, time32P/time1*100.0, stime32P, stime32P/stime1*100.0, numsolve32P, time32S/time32P);
 
   if ( COI_Error ) {
      c_log( "Errors encountered during first solve", COI_Error); }
   else if ( stacalls == 0 || solcalls == 0 ) {
      c_log( "Status or Solution routine was not during first solve called",-1 ); }
   else if ( numsolve8S != numsolve8P ) {
      c_log( "Number of solves with workpool 8 differs between static and dynamic", -1 ); }
   else if ( numsolve16S != numsolve16P ) {
      c_log( "Number of solves with workpool 16 differs between static and dynamic", -1 ); }
   else if ( numsolve32S != numsolve32P ) {
      c_log( "Number of solves with workpool 32 differs between static and dynamic", -1 ); }
 
   for ( thread=0; thread<maxthread; thread++ )
   {
      CntVect = CntVect1[thread];
      coiFree(&CntVect);
   }
   free(CntVect1);
   coiFinalize();
   c_log( "Successful Solve", OK );
}