OsiDylpSolverInterface.cpp

#ifdef COIN_USE_DYLP

#ifdef _MSC_VER

/* Turn off compiler warning about long names */
#  pragma warning(disable:4786)

/*
  MS C++ doesn't want to cope with this set of templates. Since they're just
  paranoid checks, it's easiest to simply disable them.
*/

#  define assert_same

/*
  In general, templates do not seem to be a strong point for MS C++. To get
  around this, there are a number of macros defined below.  For Sun Workshop
  and Gnu programming environments, these will expand to template functions.
  For MS, they expand to some code plus calls to STL functions. Apparently MS
  C++ doesn't have trouble with the definitions of the templates, just
  instantiations.
*/

#endif

/* Cut name lengths for readability. */

#define ODSI OsiDylpSolverInterface
#define OSI OsiSolverInterface

namespace {
  char sccsid[] = "%W%  %G%" ;
  char cvsid[] = "$Id: OsiDylpSolverInterface.cpp,v 1.22 2003/10/27 18:55:31 jpfasano Exp $" ;
}

#include <string>
#include <cassert>
#include <OsiColCut.hpp>
#include <OsiRowCut.hpp>
#include <OsiRowCutDebugger.hpp>
#include "OsiDylpSolverInterface.hpp"
#include "OsiDylpMessages.hpp"
#include "OsiDylpWarmStartBasis.hpp"

using std::string ;
using std::vector ;


extern "C"
{

#include "bonsai.h"

/*
  Dylp uses HUGE_VAL for infinity; generally this seems to resolve to IEEE
  infinity but if it doesn't, you'll need to have a look at vector.h for the
  appropriate compile-time symbol definition. A macro seems to be necessary
  here because I can't figure out how to call ODSI.getInfinity() from inside
  a static member function without explicitly passing in the value of `this'
  from the context of the call. That seems like overkill.
*/

#define DYLP_INFINITY HUGE_VAL


#ifndef ERRMSGDIR
# define ERRMSGDIR "."
#endif


extern bool dy_mpsin(const char *filename, consys_struct **consys) ;
extern void dy_initbasis(int concnt, int factor_freq, double zero_tol) ;
extern void dy_freebasis() ;

/*
  These routines work with the literal tree in DylpLib (a legacy, what can I
  say). DO NOT MIX strings handled with stralloc/strfree with strings handled
  using malloc/free.  Needless to say, std::string is another beast
  entirely.
*/

extern char *stralloc(const char *str) ;
extern void strfree(const char *str) ;

}

ioid cmdchn = IOID_NOSTRM,
     logchn = IOID_NOSTRM ;

bool cmdecho = false,
     gtxecho = false ;


int ODSI::reference_count = 0 ;
bool ODSI::basis_ready = false ;
ODSI *ODSI::dylp_owner = 0 ;





inline int ODSI::idx (int i) { return i+1 ; }

template<class T> inline T* ODSI::idx_vec (T* vec) { return vec-1 ; }

inline int ODSI::inv (int i) { return i-1 ; }

template<class T> inline T* ODSI::inv_vec (T* vec) { return vec+1 ; }

/*
  MS C++ doesn't like the template functions above. To keep it happy, wrap
  them in macros which will accomplish the same thing.
*/

#ifdef _MSC_VER
#  define IDX_VEC(zz_type_zz,zz_vec_zz) (zz_vec_zz-1)
#  define INV_VEC(zz_type_zz,zz_vec_zz) (zz_vec_zz+1)
#else
#  define IDX_VEC(zz_type_zz,zz_vec_zz) idx_vec<zz_type_zz>(zz_vec_zz)
#  define INV_VEC(zz_type_zz,zz_vec_zz) inv_vec<zz_type_zz>(zz_vec_zz)
#endif

pkvec_struct* ODSI::packed_vector (const CoinShallowPackedVector src,
                                   int dimension)

{ int n = src.getNumElements() ;
  pkvec_struct* dst = pkvec_new(n) ;

  assert(dst) ;

  if (n == 0) return (dst) ;

  packed_vector(src,dimension,dst) ;

  return (dst) ; }


void ODSI::packed_vector (const CoinShallowPackedVector src, int dimension,
                          pkvec_struct *dst)

{ int n = src.getNumElements() ;

  dst->cnt = n ;
  dst->dim = dimension ;

  if (n == 0) return ;

  const int* indices = src.getIndices() ;
  const double* elements = src.getElements() ;
  pkcoeff_struct* coeffs = dst->coeffs ;

  for (int i = 0 ; i < n ; i++)
  { coeffs[i].ndx = idx(indices[i]) ;
    coeffs[i].val = elements[i] ; }

  return ; }






/* Templates for simple copies */

template<class T> inline T* ODSI::copy (const T* src)

{ if (!src) return 0 ;
  T* dst = new T ;
  memcpy(dst, src, sizeof(T)) ;
  return dst ;
}

#ifdef _MSC_VER
#  define CLONE(zz_type_zz,zz_src_zz,zz_dst_zz) \
            if (zz_src_zz) \
            { zz_dst_zz = new zz_type_zz ; \
              memcpy(zz_dst_zz,zz_src_zz,sizeof(zz_type_zz)) ; } \
            else \
            { zz_dst_zz = NULL ; }
#else
#  define CLONE(zz_type_zz,zz_src_zz,zz_dst_zz) \
            zz_dst_zz = copy<zz_type_zz>(zz_src_zz) ;
#endif
  

template<class T> inline void ODSI::copy (const T* src, T* dst, int n)

{ if (!dst || !src || n == 0) return ;
  int size = sizeof(T) * n ;
  memcpy(dst,src,size) ;
}

#ifdef _MSC_VER
# define COPY_VEC(zz_type_zz,zz_src_zz,zz_dst_zz,zz_n_zz) \
           if (zz_dst_zz && zz_src_zz && zz_n_zz > 0 ) \
             std::copy(zz_src_zz,zz_src_zz+zz_n_zz,zz_dst_zz)
#else
# define COPY_VEC(zz_type_zz,zz_src_zz,zz_dst_zz,zz_n_zz) \
           copy<zz_type_zz>(zz_src_zz,zz_dst_zz,zz_n_zz)
#endif


template<class T> inline T* ODSI::copy (const T* src, int n)

{ if (!src || n == 0) return 0 ;
  T* dst = new T[n] ;
  copy(src, dst, n) ;
  return dst ;
}

#ifdef _MSC_VER
#  define CLONE_VEC(zz_type_zz,zz_src_zz,zz_dst_zz,zz_n_zz) \
            if (zz_src_zz && zz_n_zz > 0) \
            { zz_dst_zz = new zz_type_zz[zz_n_zz] ; \
              std::copy(zz_src_zz,zz_src_zz+zz_n_zz,zz_dst_zz) ; } \
            else \
            { zz_dst_zz = NULL ; }
#else
#  define CLONE_VEC(zz_type_zz,zz_src_zz,zz_dst_zz,zz_n_zz) \
            zz_dst_zz = copy<zz_type_zz>(zz_src_zz,zz_n_zz)
#endif


/* Specialised copy functions for complex dylp structures. */


inline void ODSI::copy_basis (const basis_struct* src, basis_struct* dst)

{ if (!src) return ;

  dst->len = src->len ;
  if (dst->el == 0)
  { throw CoinError("No basis element vector",
                    "copy_basis","OsiDylpSolverInterface") ; }
  memcpy(dst->el,src->el,idx(src->len)*sizeof(basisel_struct)) ;

# ifndef NDEBUG
  assert_same(*dst, *src, true) ;
# endif

  return ; }


inline basis_struct* ODSI::copy_basis (const basis_struct* src, int dstsze)

{ if (!src) return 0 ;
  basis_struct* dst = new basis_struct ;
  dst->el = (basisel_struct *) CALLOC(idx(dstsze),sizeof(basisel_struct)) ;
  copy_basis(src,dst) ;
  return dst ; }


lpprob_struct* ODSI::copy_lpprob (const lpprob_struct* src)

{ if (!src) return 0 ;

  int col_count = idx(src->colsze) ;
  int row_count = idx(src->rowsze) ;

  lpprob_struct* dst = NULL;
  CLONE(lpprob_struct,src,dst);

  dst->basis = copy_basis(src->basis,row_count) ;
  CLONE_VEC(flags,src->status,dst->status,col_count);
  CLONE_VEC(double,src->x,dst->x,row_count);
  CLONE_VEC(double,src->y,dst->y,row_count);
  CLONE_VEC(bool,src->actvars,dst->actvars,col_count);

# ifndef NDEBUG
  assert_same(*dst, *src, true) ;
# endif

  return dst ; }




/* Cache functions

   These properly belong in the DestructorHelpers group, but they need to be
   up here so that they are declared before the first use.
*/

inline void ODSI::destruct_col_cache ()

{ delete [] _col_x ; _col_x = 0 ;
  delete [] _col_obj ; _col_obj = 0 ;
  delete [] _col_cbar ; _col_cbar = 0 ;
}


void ODSI::destruct_row_cache ()

{ delete [] _row_lhs ; _row_lhs = 0 ;
  delete [] _row_lower ; _row_lower = 0 ;
  delete [] _row_price ; _row_price = 0 ;
  delete [] _row_range ; _row_range = 0 ;
  delete [] _row_rhs ; _row_rhs = 0 ;
  delete [] _row_sense ; _row_sense = 0 ;
  delete [] _row_upper ; _row_upper = 0 ;
}

inline void ODSI::destruct_cache ()

{ destruct_row_cache() ;
  destruct_col_cache() ;

  delete _matrix_by_row ; _matrix_by_row = 0 ;
  delete _matrix_by_col ; _matrix_by_col = 0 ;

  return ; }





inline contyp_enum ODSI::bound_to_type (double lower, double upper)

{ 

# include <cfloat>
  double inf = DBL_MAX ;

  if (upper == lower) return contypEQ ;
  bool finite_low = (lower > -inf) ;
  bool finite_up = (upper < inf) ;
  if (finite_low && finite_up) return contypRNG ;
  if (finite_low) return contypGE ;
  if (finite_up) return contypLE ;
  return contypNB ; }


inline contyp_enum ODSI::sense_to_type (char sense)

{ switch (sense)
  { case 'E': return contypEQ ;
    case 'G': return contypGE ;
    case 'L': return contypLE ;
    case 'N': return contypNB ;
    case 'R': return contypRNG ;
    default: { assert(0) ; return contypINV ; } } }


inline char ODSI::type_to_sense (contyp_enum ctypi)

{ switch (ctypi)
  { case contypNB:
    { return 'N' ; }
    case contypGE:
    { return 'G' ; }
    case contypEQ:
    { return 'E' ; }
    case contypLE:
    { return 'L' ; }
    case contypRNG:
    { return 'R' ; }
    case contypINV:
    default: 
    { assert(0) ;
      return '?' ; } } }


inline void
ODSI::gen_rowiparms (contyp_enum* ctypi, double* rhsi, double* rhslowi,
                     char sensei, double rhsini, double rangei)

{ *ctypi = sense_to_type(sensei) ;
  switch (*ctypi)
  { case contypEQ:
    { *rhslowi = 0 ;
      *rhsi = rhsini ;
      break ; }
    case contypLE:
    { *rhslowi = 0 ;
      *rhsi = rhsini ;
      break ; }
    case contypGE:
    { *rhslowi = 0 ;
      *rhsi = rhsini ;
      break ; }
    case contypRNG:
    { *rhslowi = rhsini-rangei ;
      *rhsi = rhsini ;
      break ; }
    case contypNB:
    { *rhslowi = 0 ;
      *rhsi = 0 ;
      break ; }
    default:
    { assert(false) ; } }

  return ; }


inline void
ODSI::gen_rowiparms (contyp_enum* ctypi, double* rhsi, double* rhslowi,
                     double rowlbi, double rowubi)

{ *ctypi = bound_to_type(rowlbi,rowubi) ;
  switch (*ctypi)
  { case contypEQ:
    case contypLE:
    { *rhslowi = 0 ;
      *rhsi = rowubi ;
      break ; }
    case contypGE:
    { *rhslowi = 0 ;
      *rhsi = rowlbi ;
      break ; }
    case contypRNG:
    { *rhslowi = rowlbi ;
      *rhsi = rowubi ;
      break ; }
    case contypNB:
    { *rhslowi = 0 ;
      *rhsi = 0 ;
      break ; }
    default:
    { assert(false) ; } }

  return ; }


void ODSI::add_col (const CoinPackedVectorBase& coin_colj,
                    vartyp_enum vtypj, double vlbj, double vubj, double objj)

{ pkvec_struct *pk_colj = packed_vector(coin_colj,getNumRows()) ;

/*
  Do we need a consys?
*/
  if (!consys) construct_consys(0,0) ;
/*
  Add the column.
*/
  bool r = consys_addcol_pk(consys,vtypj,pk_colj,objj,vlbj,vubj) ;
  pkvec_free(pk_colj) ;
  assert(r) ;
/*
  After adding a column, the best we can do is a warm start.
*/
  resolveOptions->forcewarm = true ;

  destruct_cache() ; }


void ODSI::add_row (const CoinPackedVectorBase &coin_rowi, char clazzi,
                    contyp_enum ctypi, double rhsi, double rhslowi)

{ pkvec_struct *pk_rowi = packed_vector(coin_rowi,getNumCols()) ;

/*
  Do we need a consys?
*/
  if (!consys) construct_consys(0,0) ;
/*
  Add the row.
*/
  bool r = consys_addrow_pk(consys,clazzi,ctypi,pk_rowi,rhsi,rhslowi,0,0) ;
  pkvec_free(pk_rowi) ;
  assert(r) ;
/*
  After adding a constraint, the best we can do is a warm start.
*/
  resolveOptions->forcewarm = true ;

  destruct_cache() ; }


void ODSI::worst_case_primal ()

{ 

/*
  No columns, no solution. A non-zero value guarantees a consys structure.
*/
  int n = getNumCols() ;
  if (n == 0) return ;
  assert(consys && consys->obj && consys->vub && consys->vlb) ;

  double *obj = consys->obj ;
  double *vlb = consys->vlb ;
  double *vub = consys->vub ;
/*
  We have columns. Allocate a cached primal solution vector, if it doesn't
  already exist.
*/
  if (!_col_x) _col_x = new double[n] ;
/*
  Walk the objective vector, taking values for variables from the appropriate
  bounds vector. The objective attached to consys is always a minimisation
  objective.
*/
  for (int j = 1 ; j <= n ; j++)
  { if (obj[j] > 0)
    { _col_x[inv(j)] = vub[j] ; }
    else
    { _col_x[inv(j)] = vlb[j] ; } }
  
  return ; }

void ODSI::calc_objval ()

{ int n = getNumCols() ;

/*
  The easy case --- if we have no variables, the objective must be 0.
*/
  if (n == 0)
  { _objval = 0 ;
    return ; }
/*
  We have variables. Grab the primal solution and objective function and
  calculate their dot product.
*/
  const double *obj = getObjCoefficients() ;
  const double *x = getColSolution() ;
  _objval = 0.0 ;
  for (int j = 0 ; j < n ; j++) _objval += x[j]*obj[j] ;
  setcleanzero(_objval,tolerances->cost) ;

  return ; }






/* OsiDylpSolverInterface constructors and related subroutines. */



inline void ODSI::construct_lpprob ()

{ lpprob = new lpprob_struct ;
  memset(lpprob, 0, sizeof(lpprob_struct)) ;
  setflg(lpprob->ctlopts,lpctlNOFREE) ;
  lpprob->phase = dyINV ;
  lpprob->consys = consys ;
  lpprob->rowsze = consys->rowsze ;
  lpprob->colsze = consys->colsze ;
  
  return ; }


inline void ODSI::construct_options ()

{ 
/*
  Acquire the default options and tolerances from dylp. Set the OSI options
  and tolerances to match. Turn dylp's printing options down to `catatonic'.
*/
  delete initialSolveOptions ;
  initialSolveOptions = new lpopts_struct ;
  delete tolerances ;
  tolerances = new lptols_struct ;
  dy_defaults(&initialSolveOptions,&tolerances) ;
  delete resolveOptions ;
  CLONE(lpopts_struct,initialSolveOptions,resolveOptions);
  dy_setprintopts(0,initialSolveOptions) ;
  dy_setprintopts(0,resolveOptions) ;

  setIntParam(OsiMaxNumIteration,3*initialSolveOptions->iterlim) ;
  setIntParam(OsiMaxNumIterationHotStart,3*initialSolveOptions->iterlim) ;
  setDblParam(OsiDualTolerance,tolerances->dfeas_scale*tolerances->cost) ;
  setDblParam(OsiPrimalTolerance,tolerances->pfeas_scale*tolerances->zero) ;
/*
  Differentiate options for initial solution vs. reoptimising.
*/
  initialSolveOptions->forcecold = true ;
  initialSolveOptions->fullsys = true ;
  resolveOptions->forcecold = false ;
  resolveOptions->fullsys = false ;
/*
  Acquire and clear a statistics structure.
*/
  delete statistics ;
  statistics = new lpstats_struct ;
  memset(statistics,0,sizeof(lpstats_struct)) ;
}

void ODSI::construct_consys (int cols, int rows)

{
/*
  Set parameters to specify the appropriate options and attached vectors.
  
  parts specifies (in order) an objective, variable upper bounds, variable
  lower bounds, constraint rhs, range constraint rhs, variable type, and
  constraint type.

  options specifies that consys should issue a warning if requested to attach
  a vector that's already attached.
*/
  flags parts = CONSYS_OBJ | 
                CONSYS_VUB | CONSYS_VLB | CONSYS_VTYP |
                CONSYS_RHS | CONSYS_RHSLOW | CONSYS_CTYP ;
  flags opts = CONSYS_WRNATT ;
/*
  The first parameter to consys is the constraint system name, which is
  never conveniently available in the OSI frame.
*/
  consys = consys_create(0,parts,opts,rows,cols) ;
  assert(consys) ;

  return ; }


void ODSI::dylp_ioinit ()

{ if (reference_count > 1) return ;

  std::string errfile = std::string(ERRMSGDIR)+std::string("/bonsaierrs.txt") ;
# ifndef NDEBUG
  errinit(const_cast<char *>(errfile.c_str()),0,true) ;
# else
  errinit(const_cast<char *>(errfile.c_str()),0,false) ;
# endif
  bool r1 = ioinit() ;
  assert(r1) ;
}

void ODSI::load_problem (const CoinMpsIO &mps)

{ int n = mps.getNumCols() ;
  int m = mps.getNumRows() ;
/*
  Unpack the operands from the CoinMpsIO object.
*/
  const CoinPackedMatrix matrix = *mps.getMatrixByCol() ;
  const double* col_lower = mps.getColLower() ;
  const double* col_upper = mps.getColUpper() ;
  const double* obj = mps.getObjCoefficients() ;
  const char *sense = mps.getRowSense() ;
  const double* rhsin = mps.getRightHandSide() ;
  const double* range = mps.getRowRange() ;

  double *rhs = new double[m] ;
  double *rhslow = new double[m] ; 
  contyp_enum *ctyp = new contyp_enum[m] ;

  gen_rowparms(m,rhs,rhslow,ctyp,sense,rhsin,range) ;
/*
  Free the existing problem structures, preserving options and tolerances
  and resetting statistics.
*/
  destruct_problem(true) ;
/*
  Create an empty consys_struct. Load the constraint system name, objective
  name, and objective offset.
*/
  construct_consys(n,m) ;
  setStrParam(OsiProbName,mps.getProblemName()) ;
  if (consys->objnme) strfree(consys->objnme) ;
  consys->objnme = stralloc(mps.getObjectiveName()) ;
  setDblParam(OsiObjOffset,mps.objectiveOffset()) ;
/*
  First loop: Insert empty constraints into the new constraint system.
*/
  pkvec_struct* rowi = pkvec_new(0) ;
  assert(rowi) ;

  for (int i = 0 ; i < m ; i++)
  { rowi->nme = const_cast<char *>(mps.rowName(i)) ;
    bool r = consys_addrow_pk(consys,'a',ctyp[i],rowi,rhs[i],rhslow[i],0,0) ;
    assert(r) ; }
  
  if (rowi) pkvec_free(rowi) ;
  delete[] rhs ;
  delete[] rhslow ;
  delete[] ctyp ;
/*
  Take a moment and set up a vector of variable types. It's a bit of a pain
  that CoinMpsIO doesn't distinguish binary from general integer, but we can
  do it here by checking bounds. The test is drawn from OSI::isBinary and
  will include variables fixed at 0 or 1.
*/
  const char *const intvars = mps.integerColumns() ;
  vartyp_enum *const vtyp = new vartyp_enum[n] ;
  if (intvars)
  { for (int j = 0 ; j < n ; j++)
    { if (intvars[j])
      { if ((col_lower[j] == 0.0 || col_lower[j] == 1.0) &&
            (col_upper[j] == 0.0 || col_upper[j] == 1.0))
        { vtyp[j] = vartypBIN ; }
        else
        { vtyp[j] = vartypINT ; } }
      else
      { vtyp[j] = vartypCON ; } } }
  else
  { for (int j = 0 ; j < n ; j++) vtyp[j] = vartypCON ; }
/*
  Second loop. Insert the coefficients by column. If we need to create a
  column-ordered copy, take advantage and cache it.  The size of colj is a
  gross overestimate, but it saves the trouble of constantly reallocating it.
*/
  const CoinPackedMatrix *matrix2 ;
  if (matrix.isColOrdered())
  { matrix2 = &matrix ; }
  else
  { _matrix_by_col = new CoinPackedMatrix ;
    _matrix_by_col->reverseOrderedCopyOf(matrix) ;
    matrix2 = _matrix_by_col ; }
  
  pkvec_struct* colj = pkvec_new(n) ;

  for (int j = 0 ; j < n ; j++)
  { const CoinShallowPackedVector coin_col = matrix2->getVector(j) ;
    packed_vector(coin_col,n,colj) ;
    colj->nme = const_cast<char *>(mps.columnName(j)) ;
    bool r = consys_addcol_pk(consys,vtyp[j],colj,obj[j],
                              col_lower[j],col_upper[j]) ;
    assert(r) ; }

  pkvec_free(colj) ;
  delete[] vtyp ;
  assert(matrix2->isEquivalent(*getMatrixByCol())) ;

/*
  Finish up. Establish a primal solution and an objective. The primal solution
  is `worst-case', in the sense that it's constructed by setting each primal
  variable to the bound that maximises the objective.
*/
  worst_case_primal() ;
  calc_objval() ;

  return ; }

void ODSI::load_problem (const CoinPackedMatrix& matrix,
    const double* col_lower, const double* col_upper, const double* obj,
    const contyp_enum *ctyp, const double* rhs, const double* rhslow)

{ 
/*
  Free the existing problem structures, preserving options and tolerances
  and resetting statistics.
*/
  destruct_problem(true) ;
/*
  Create an empty consys_struct.
*/
  int m = matrix.getNumCols() ;
  int n = matrix.getNumRows() ;

  construct_consys(n,m) ;
/*
  First loop: Insert empty constraints into the new constraint system.
*/
  pkvec_struct* rowi = pkvec_new(0) ;
  assert(rowi) ;

  for (int i = 0 ; i < n ; i++)
  { rowi->nme = 0 ;
    bool r = consys_addrow_pk(consys,'a',ctyp[i],rowi,rhs[i],rhslow[i],0,0) ;
    assert(r) ;
  }

  if (rowi) pkvec_free(rowi) ;
/*
  Second loop. Insert the coefficients by column. If we need to create a
  column-ordered copy, take advantage and cache it.  The size of colj is a
  gross overestimate, but it saves the trouble of constantly reallocating it.
*/
  const CoinPackedMatrix *matrix2 ;
  if (matrix.isColOrdered())
  { matrix2 = &matrix ; }
  else
  { _matrix_by_col = new CoinPackedMatrix ;
    _matrix_by_col->reverseOrderedCopyOf(matrix) ;
    matrix2 = _matrix_by_col ; }
  
  pkvec_struct* colj = pkvec_new(n) ;

  for (int j = 0 ; j < m ; j++)
  { const CoinShallowPackedVector coin_col = matrix2->getVector(j) ;
    packed_vector(coin_col,n,colj) ;
    double objj = obj?obj[j]:0 ;
    double vlbj = col_lower?col_lower[j]:0 ;
    double vubj = col_upper?col_upper[j]:DYLP_INFINITY ;
    colj->nme = 0 ;
    bool r = consys_addcol_pk(consys,vartypCON,colj,objj,vlbj,vubj) ;
    assert(r) ; }

  pkvec_free(colj) ;
  assert(matrix2->isEquivalent(*getMatrixByCol())) ;

/*
  Finish up. Establish a primal solution and an objective. The primal solution
  is `worst-case', in the sense that it's constructed by setting each primal
  variable to the bound that maximises the objective.
*/
  worst_case_primal() ;
  calc_objval() ;

  return ; }


void ODSI::load_problem (const int colcnt, const int rowcnt,
           const int *start, const int *index, const double *value,
           const double* col_lower, const double* col_upper, const double* obj,
           const contyp_enum *ctyp, const double* rhs, const double* rhslow)

{ 
/*
  Free the existing problem structures, preserving options and tolerances
  and resetting statistics.
*/
  destruct_problem(true) ;
/*
  Create an empty consys_struct. Request the standard attached vectors:
  objective, variable upper & lower bounds, variable and constraint types,
  and constraint right-hand-side and range (rhslow).
*/
  construct_consys(rowcnt,colcnt) ;
/*
  First loop: Insert empty constraints into the new constraint system.
*/
  pkvec_struct* rowi = pkvec_new(0) ;
  assert(rowi) ;

  for (int i = 0 ; i < rowcnt ; i++)
  { rowi->nme = 0 ;
    bool r = consys_addrow_pk(consys,'a',ctyp[i],rowi,rhs[i],rhslow[i],0,0) ;
    assert(r) ;
  }

  if (rowi) pkvec_free(rowi) ;
/*
  Second loop. Insert the coefficients by column. The size of colj is a gross
  overestimate, but it saves the trouble of constantly reallocating it.
*/
  pkvec_struct *colj = pkvec_new(rowcnt) ;
  assert(colj) ;
  colj->dim = rowcnt ;
  pkcoeff_struct *coeffs = colj->coeffs ;

  for (int j = 0 ; j < colcnt ; j++)
  { int startj = start[j] ;
    int lenj = start[j+1]-startj ;

    for (int ndx = 0 ; ndx < lenj ; ndx++)
    { coeffs[ndx].ndx = idx(index[startj+ndx]) ;
      coeffs[ndx].val = value[startj+ndx] ; }
    colj->cnt = lenj ;
  
    double objj = obj?obj[j]:0 ;
    double vlbj = col_lower?col_lower[j]:0 ;
    double vubj = col_upper?col_upper[j]:DYLP_INFINITY ;
    colj->nme = 0 ;
    bool r = consys_addcol_pk(consys,vartypCON,colj,objj,vlbj,vubj) ;
    assert(r) ;
  }

  if (colj) pkvec_free(colj) ;
/*
  Finish up. Establish a primal solution and an objective.
*/
  worst_case_primal() ;
  calc_objval() ;
  
  return ; }



void ODSI::gen_rowparms (int rowcnt,
                         double *rhs, double *rhslow, contyp_enum *ctyp,
                         const double *rowlb, const double *rowub)

{ double inf = DYLP_INFINITY ;
  double rowlbi,rowubi ;

  for (int i = 0 ; i < rowcnt ; i++)
  { rowlbi = rowlb?rowlb[i]:-inf ;
    rowubi = rowub?rowub[i]:inf ;
    ctyp[i] = bound_to_type(rowlbi,rowubi) ;

    switch (ctyp[i])
    { case contypEQ:
      case contypLE:
      { rhs[i] = rowubi ; 
        rhslow[i] = 0 ;
        break ; }
      case contypGE:
      { rhs[i] = rowlbi ;
        rhslow[i] = 0 ;
        break ; }
      case contypRNG:
      { rhs[i] = rowubi ; 
        rhslow[i] = rowlbi ;
        break ; }
      case contypNB:
      { rhs[i] = inf ;
        rhslow[i] = -inf ;
        break ; }
      default:
      { assert(0) ; } } }

  return ; }


void ODSI::gen_rowparms (int rowcnt,
                         double *rhs, double *rhslow, contyp_enum *ctyp,
                 const char *sense, const double *rhsin, const double *range)

{ double rhsi,rangei ;
  char sensei ;

  for (int i = 0 ; i < rowcnt ; i++)
  { rhsi = rhsin?rhsin[i]:0 ;
    rangei = range?range[i]:0 ;
    sensei = sense?sense[i]:'G' ;
    
    switch (sensei)
    { case 'E':
      { rhs[i] = rhsi ;
        rhslow[i] = 0 ;
        ctyp[i] = contypEQ ;
        break ; }
      case 'L':
      { rhs[i] = rhsi ;
        rhslow[i] = 0 ;
        ctyp[i] = contypLE ;
        break ; }
      case 'G':
      { rhs[i] = rhsi ;
        rhslow[i] = 0 ;
        ctyp[i] = contypGE ;
        break ; }
      case 'R':
      { rhs[i] = rhsi ;
        rhslow[i] = rhsi-rangei ;
        ctyp[i] = contypRNG ;
        break ; }
      case 'N':
      { rhs[i] = 0 ;
        rhslow[i] = 0 ;
        ctyp[i] = contypNB ;
        break ; }
      default:
      { assert(0) ; } } }

  return ; }





ODSI::OsiDylpSolverInterface ()

  : initialSolveOptions(0),
    resolveOptions(0),
    tolerances(0),
    consys(0),
    lpprob(0),
    statistics(0),

    local_logchn(IOID_NOSTRM),
    initial_gtxecho(false),
    resolve_gtxecho(false),
    lp_retval(lpINV),
    obj_sense(1.0),

    solvername("dylp"),
    mps_debug(0),
    hotstart_fallback(0),
    activeBasis(0),
    activeIsModified(false),

    _objval(0),
    _col_x(0),
    _col_obj(0),
    _col_cbar(0),
    _row_lhs(0),
    _row_lower(0),
    _row_price(0),
    _row_range(0),
    _row_rhs(0),
    _row_sense(0),
    _row_upper(0),
    _matrix_by_row(0),
    _matrix_by_col(0)

{
/*
  Replace the OSI default messages with ODSI messages.
*/
  setOsiDylpMessages(CoinMessages::us_en) ;
/*
  Clear the hint info_ array.
*/
  for (int i = 0 ; i < OsiLastHintParam ; i++) info_[i] = 0 ;
/*
  Acquire the default options and tolerances.
*/
  construct_options() ;
/*
  Increment the ODSI existence count. If this is the first interface in
  existence, initialise the i/o package.
*/
  reference_count++ ;

  if (reference_count == 1)
  { dylp_ioinit() ;
    CoinRelFltEq eq ;
    assert(eq(DYLP_INFINITY, DYLP_INFINITY)) ; }
  
  return ; }


ODSI::OsiDylpSolverInterface (const OsiDylpSolverInterface& src)

  : OsiSolverInterface(src),
    initial_gtxecho(src.initial_gtxecho),
    resolve_gtxecho(src.resolve_gtxecho),
    lp_retval(src.lp_retval),
    obj_sense(src.obj_sense),

    solvername(src.solvername),
    mps_debug(src.mps_debug),
    hotstart_fallback(0),       // do not copy hot start information
    activeBasis(0),
    activeIsModified(false),
  
    _objval(src._objval),
    _matrix_by_row(0),
    _matrix_by_col(0)

{ if (src.consys)
  { bool r = consys_dupsys(src.consys,&consys,src.consys->parts) ;
    assert(r) ; }
  else
  { consys = 0 ; }
  if (src.lpprob)
  { lpprob = copy_lpprob(src.lpprob) ;
    lpprob->consys = consys ; }
  else
  { lpprob = 0 ; }

  CLONE(lpopts_struct,src.initialSolveOptions,initialSolveOptions) ;
  CLONE(lpopts_struct,src.resolveOptions,resolveOptions) ;
  CLONE(lpstats_struct,src.statistics,statistics) ;
  CLONE(lptols_struct,src.tolerances,tolerances) ;

  if (src.local_logchn != IOID_NOSTRM)
  { local_logchn = openfile(idtopath(src.local_logchn),
                            const_cast<char *>("w")) ; }
  else
  { local_logchn = IOID_NOSTRM ; }
  assert(local_logchn == src.local_logchn) ;

  if (src.activeBasis)
  { activeBasis = src.activeBasis->clone() ;
    activeIsModified = src.activeIsModified ; }

  int col_count = src.getNumCols() ;
  int row_count = src.getNumRows() ;

  CLONE_VEC(double,src._col_x,_col_x,col_count) ;
  CLONE_VEC(double,src._col_obj,_col_obj,col_count) ;
  CLONE_VEC(double,src._col_cbar,_col_cbar,col_count) ;
  CLONE_VEC(double,src._row_lhs,_row_lhs,row_count) ;
  CLONE_VEC(double,src._row_lower,_row_lower,row_count) ;
  CLONE_VEC(double,src._row_price,_row_price,row_count) ;
  CLONE_VEC(double,src._row_range,_row_range,row_count) ;
  CLONE_VEC(double,src._row_rhs,_row_rhs,row_count) ;
  CLONE_VEC(char,src._row_sense,_row_sense,row_count) ;
  CLONE_VEC(double,src._row_upper,_row_upper,row_count) ;

  if (src._matrix_by_row) 
    _matrix_by_row = new CoinPackedMatrix(*src._matrix_by_row) ;
  if (src._matrix_by_col) 
    _matrix_by_col = new CoinPackedMatrix(*src._matrix_by_col) ;

  COPY_VEC(void *,&src.info_[0],&info_[0],OsiLastHintParam) ;

  reference_count++ ;

# ifndef NDEBUG
  assert_same(*this, src, true) ;
# endif

}


inline OsiSolverInterface* ODSI::clone (bool copyData) const

{ if (copyData)
  { return new OsiDylpSolverInterface(*this) ; }
  else
  { return new OsiDylpSolverInterface() ; }
}


OsiDylpSolverInterface &ODSI::operator= (const OsiDylpSolverInterface &rhs)

{ if (this != &rhs)
  { destruct_problem(false) ;
  
    OSI::operator=(rhs) ;

    if (rhs.consys)
    { bool r = consys_dupsys(rhs.consys,&consys,rhs.consys->parts) ;
      assert(r) ; }
    else
    { consys = 0 ; }
    if (rhs.lpprob)
    { lpprob = copy_lpprob(rhs.lpprob) ;
      lpprob->consys = consys ; }
    else
    { lpprob = 0 ; }

    CLONE(lpopts_struct,rhs.initialSolveOptions,initialSolveOptions) ;
    CLONE(lpopts_struct,rhs.resolveOptions,resolveOptions) ;
    CLONE(lpstats_struct,rhs.statistics,statistics) ;
    CLONE(lptols_struct,rhs.tolerances,tolerances) ;

    if (rhs.local_logchn != IOID_NOSTRM)
    { local_logchn = openfile(idtopath(rhs.local_logchn),
                              const_cast<char *>("w")) ; }
    else
    { local_logchn = IOID_NOSTRM ; }
    assert(local_logchn == rhs.local_logchn) ;

    if (rhs.activeBasis)
    { activeBasis = rhs.activeBasis->clone() ;
      activeIsModified = rhs.activeIsModified ; }

    int col_count = rhs.getNumCols() ;
    int row_count = rhs.getNumRows() ;

    CLONE_VEC(double,rhs._col_x,_col_x,col_count) ;
    CLONE_VEC(double,rhs._col_obj,_col_obj,col_count) ;
    CLONE_VEC(double,rhs._col_cbar,_col_cbar,col_count) ;
    CLONE_VEC(double,rhs._row_lhs,_row_lhs,row_count) ;
    CLONE_VEC(double,rhs._row_lower,_row_lower,row_count) ;
    CLONE_VEC(double,rhs._row_price,_row_price,row_count) ;
    CLONE_VEC(double,rhs._row_range,_row_range,row_count) ;
    CLONE_VEC(double,rhs._row_rhs,_row_rhs,row_count) ;
    CLONE_VEC(char,rhs._row_sense,_row_sense,row_count) ;
    CLONE_VEC(double,rhs._row_upper,_row_upper,row_count) ;

    if (rhs._matrix_by_row) 
      _matrix_by_row = new CoinPackedMatrix(*rhs._matrix_by_row) ;
    if (rhs._matrix_by_col) 
      _matrix_by_col = new CoinPackedMatrix(*rhs._matrix_by_col) ;

    COPY_VEC(void *,&rhs.info_[0],&info_[0],OsiLastHintParam) ;

    reference_count++ ;

#   ifndef NDEBUG
    assert_same(*this, rhs, true) ;
#   endif
  }

  return (*this) ; }





/* OsiDylpSolverInterface destructor and related subroutines. */


void ODSI::destruct_problem (bool preserve_interface)

{ 
/*
  If this object claims ownership of the solver, it should possess an lpprob
  structure created when the solver was called.
*/
  assert((dylp_owner != this) || (dylp_owner == this && lpprob)) ;

  if (dylp_owner == this) detach_dylp() ;

  if (lpprob)
  { assert(lpprob->consys == consys) ;
    consys_free(consys) ;
    consys = 0 ;
    dy_freesoln(lpprob) ;
    delete lpprob ;
    lpprob = 0 ; }
  else
  if (consys)
  { consys_free(consys) ;
    consys = 0 ; }
 
  if (hotstart_fallback)
  { delete hotstart_fallback ;
    hotstart_fallback = 0 ; }
  if (activeBasis)
  { delete activeBasis ;
    activeBasis = 0 ;
    activeIsModified = false ; }

  destruct_cache() ;

  if (preserve_interface == false)
  { if (initialSolveOptions)
    { delete initialSolveOptions ;
      initialSolveOptions = 0 ; }
    if (resolveOptions)
    { delete resolveOptions ;
      resolveOptions = 0 ; }
    if (tolerances)
    { delete tolerances ;
      tolerances = 0 ; } }
  
  if (statistics)
  { if (preserve_interface == true)
    { memset(statistics,0,sizeof(lpstats_struct)) ; }
    else
    { delete statistics ;
      statistics = 0 ; } }
  
  return ; }


void ODSI::detach_dylp ()

{ assert(dylp_owner == this && lpprob && lpprob->consys) ;
  assert(flgon(lpprob->ctlopts,lpctlDYVALID)) ;
  flags save_flags = getflg(lpprob->ctlopts,lpctlNOFREE|lpctlONLYFREE) ;
  clrflg(lpprob->ctlopts,lpctlNOFREE) ;
  setflg(lpprob->ctlopts,lpctlONLYFREE) ;
  lpprob->phase = dyDONE ;
/*
  Either of the initialSolve or resolve option blocks are ok here; dylp does
  not look.
*/
  dylp(lpprob,initialSolveOptions,tolerances,statistics) ;
  clrflg(lpprob->ctlopts,lpctlONLYFREE) ;
  setflg(lpprob->ctlopts,save_flags) ;
  dylp_owner = 0 ; }


ODSI::~OsiDylpSolverInterface ()

{
/*
  Destroy the problem-specific structures used by dylp, and close the log
  file, if open.
*/
  destruct_problem(false) ;
  if (local_logchn != IOID_INV && local_logchn != IOID_NOSTRM)
    (void) closefile(local_logchn) ;

  reference_count-- ;
  if (reference_count == 0)
  { if (basis_ready == true)
    { dy_freebasis() ;
      basis_ready = false ; }
    ioterm() ;
    errterm() ; }
  
  return ; }

void ODSI::reset ()

{ 
/*
  Destroy the problem-specific structures used by dylp, and close the log
  file, if open.
*/
  destruct_problem(false) ;
  if (local_logchn != IOID_INV && local_logchn != IOID_NOSTRM)
  { (void) closefile(local_logchn) ;
    local_logchn = IOID_NOSTRM ; }
/*
  That takes care of cleaning out the ODSI object. Call setInitialData to
  reset the underlying OSI object.
*/
  setInitialData() ;
/*
  Now the equivalent for ODSI: ensure default values for various fields in the
  object. Regrettably, messages_ is not a pointer, so we can't preserve the
  CoinMessages structure over the reset.
*/
  initial_gtxecho = false ;
  resolve_gtxecho = false ;
  lp_retval = lpINV ;
  setObjSense(1.0) ;
  mps_debug = false ;
  construct_options() ;
  setOsiDylpMessages(CoinMessages::us_en) ;
  for (int i = 0 ; i < OsiLastHintParam ; i++) info_[i] = 0 ;

  return ; }




inline void ODSI::setContinuous (int j)

{ if (!consys->vtyp)
  { bool r = consys_attach(consys,CONSYS_VTYP,sizeof(vartyp_enum),
                           reinterpret_cast<void **>(&consys->vtyp)) ;
    assert(r) ; }

  consys->vtyp[idx(j)] = vartypCON ; }


inline void ODSI::setContinuous (const int* indices, int len)

{ if (!consys->vtyp)
  { bool r = consys_attach(consys,CONSYS_VTYP,sizeof(vartyp_enum),
                           reinterpret_cast<void **>(&consys->vtyp)) ;
    assert(r) ; }

  for (int i = 0 ; i < len ; i++)
    consys->vtyp[idx(indices[i])] = vartypCON ; }


inline void ODSI::setInteger (int j)

{ if (!consys->vtyp)
  { bool r = consys_attach(consys,CONSYS_VTYP,sizeof(vartyp_enum),
                           reinterpret_cast<void **>(&consys->vtyp)) ;
    assert(r) ; }

  if (getColLower()[j] == 0.0 && getColUpper()[j] == 1.0)
    consys->vtyp[idx(j)] = vartypBIN ;
  else
    consys->vtyp[idx(j)] = vartypINT ; }


inline void ODSI::setInteger (const int* indices, int len)

{ if (!consys->vtyp)
  { bool r = consys_attach(consys,CONSYS_VTYP,sizeof(vartyp_enum),
                           reinterpret_cast<void **>(&consys->vtyp)) ;
    assert(r) ; }

  for (int i = 0 ; i < len ; i++) setInteger(indices[i]) ; }


inline void ODSI::setColLower (int i, double val)

{ if (!consys->vlb)
  { bool r = consys_attach(consys,CONSYS_VLB,sizeof(double),
                           reinterpret_cast<void **>(&consys->vlb)) ;
    assert(r) ; }

  consys->vlb[idx(i)] = val ;
  if (lpprob) setflg(lpprob->ctlopts,lpctlLBNDCHG) ;

  if (isInteger(i))
  { if (floor(val) != val)
      setContinuous(i) ;
    else
    if (isBinary(i) && !(val == 0.0 || val == 1.0))
      setInteger(i) ; } }

inline void ODSI::setColUpper (int i, double val)

{ if (!consys->vub)
  { bool r = consys_attach(consys,CONSYS_VUB,sizeof(double),
                           reinterpret_cast<void **>(&consys->vub)) ;
    assert(r) ; }

  consys->vub[idx(i)] = val ;
  if (lpprob) setflg(lpprob->ctlopts,lpctlUBNDCHG) ;

  if (isInteger(i))
  { if (floor(val) != val)
      setContinuous(i) ;
    else
    if (isBinary(i) && !(val == 0.0 || val == 1.0))
      setInteger(i) ; } }


void ODSI::setRowType (int i, char sense, double rhs, double range)

{ int k = idx(i) ;

  gen_rowiparms(&consys->ctyp[k],&consys->rhs[k],&consys->rhslow[k],
                sense,rhs,range) ;
  if (resolveOptions) resolveOptions->forcewarm = true ;

  destruct_row_cache() ; }


void ODSI::setRowUpper (int i, double val)

{ int k = idx(i) ;
  double clbi ;

  switch (consys->ctyp[k])
  { case contypEQ:
    case contypGE:
    { clbi = consys->rhs[k] ;
      break ; }
    case contypRNG:
    { clbi = consys->rhslow[k] ;
      break ; }
    case contypLE:
    case contypNB:
    { clbi = -DYLP_INFINITY ;
      break ; }
    default:
    { assert(false) ; } }
  
  gen_rowiparms(&consys->ctyp[k],&consys->rhs[k],&consys->rhslow[k],
                clbi,val) ;
  if (lpprob) setflg(lpprob->ctlopts,lpctlRHSCHG) ;

  destruct_row_cache() ; }


void ODSI::setRowLower (int i, double val)

{ int k = idx(i) ;
  double cubi ;
  contyp_enum ctypi = consys->ctyp[k] ;

  if (ctypi == contypGE || ctypi == contypNB)
    cubi = DYLP_INFINITY ;
  else
    cubi = consys->rhs[k] ;

  gen_rowiparms(&consys->ctyp[k],&consys->rhs[k],&consys->rhslow[k],
                val,cubi) ;
  if (lpprob) setflg(lpprob->ctlopts,lpctlRHSCHG) ;

  destruct_row_cache() ; }


inline void ODSI::addRow (const CoinPackedVectorBase &row,
                          double rlb, double rub)

{ contyp_enum ctypi ;
  double rhsi,rhslowi ;

  gen_rowiparms(&ctypi,&rhsi,&rhslowi,rlb,rub) ;
  add_row(row,'a',ctypi,rhsi,rhslowi) ;

  return ; }


inline void ODSI::addRow (const CoinPackedVectorBase &coin_row, 
                          char sense, double rhs, double range)

{ contyp_enum ctypi ;
  double rhsi,rhslowi ;

  gen_rowiparms(&ctypi,&rhsi,&rhslowi,sense,rhs,range) ;
  add_row(coin_row,'a',ctypi,rhsi,rhslowi) ;

  return ; }


void ODSI::applyRowCut (const OsiRowCut &cut)

{ contyp_enum ctypi ;
  double rhsi,rhslowi ;

  gen_rowiparms(&ctypi,&rhsi,&rhslowi,cut.lb(),cut.ub()) ;
  add_row(cut.row(),'c',ctypi,rhsi,rhslowi) ;
  
  return ; }


void ODSI::deleteRows (int count, const int* rows)

{ if (count <= 0) return ;
/*
  Sort the indices and do the deletions, one by one.
*/
  vector<int> lclrows = vector<int>(&rows[0],&rows[count]) ;

  if (count > 1) std::sort(lclrows.begin(),lclrows.end()) ;

  for (int k = count-1 ; k >= 0 ; k--)
  { int i = idx(lclrows[k]) ;
    bool r = consys_delrow_stable(consys,i) ;
    assert(r) ; }
/*
  Now, see if there's an active basis. If so, check that all the constraints
  to be deleted are slack. If they are, we can delete them from activeBasis
  and still guarantee a valid basis. If not, throw away activeBasis.
*/
  if (activeBasis)
  { bool allslack = true ;
    OsiDylpWarmStartBasis *odwsb =
      dynamic_cast<OsiDylpWarmStartBasis *>(activeBasis) ;
    for (int k = count-1 ; k >= 0 ; k--)
    { int i = lclrows[k] ;
      if (odwsb->getArtifStatus(i) != CoinWarmStartBasis::basic)
      { allslack = false ;
        break ; } }
    if (allslack == true)
    { odwsb->deleteRows(count,rows) ;
      activeIsModified = true ;
      resolveOptions->forcewarm = true ; }
    else
    { delete activeBasis ;
      activeBasis = 0 ;
      activeIsModified = false ;
      resolveOptions->forcecold = true ; } }

  destruct_cache() ; }



void ODSI::setObjCoeff (int j, double objj)

{ if (j >= getNumCols()) return ;

  consys->obj[idx(j)] = getObjSense()*objj ;
  if (_col_obj) _col_obj[j] = objj ;
  if (lpprob) setflg(lpprob->ctlopts,lpctlOBJCHG) ;
  calc_objval() ; }
  

void ODSI::setObjSense (double val)
/*
  The `natural' action of OSI (and dylp) is minimisation. Maximisation is
  accomplished as min -cx. So using -1 for maximisation is more natural
  than you'd think at first glance.
*/

{ int n = getNumCols() ;

  if (n > 0 && val != obj_sense)
  { double *tmpobj = INV_VEC(double,consys->obj) ;
    std::transform(tmpobj,tmpobj+n,tmpobj,std::negate<double>()) ;
    if (lpprob) setflg(lpprob->ctlopts,lpctlOBJCHG) ; }
  
  obj_sense = val ;
  
  return ; }


inline void ODSI::addCol (const CoinPackedVectorBase& coin_col,
                          double vlb, double vub, double obj)

{ add_col(coin_col,vartypCON,vlb,vub,obj) ; }


void ODSI::applyColCut (const OsiColCut &cut)

{ const double* old_lower = getColLower() ;
  const double* old_upper = getColUpper() ;
  const CoinPackedVector& new_lower = cut.lbs() ;
  const CoinPackedVector& new_upper = cut.ubs() ;

  int i ;
  int n1 = new_lower.getNumElements() ;
  int n2 = new_upper.getNumElements() ;

  for (i = 0 ; i < n1 ; i++)
  { int j = new_lower.getIndices()[i] ;
    double x = new_lower.getElements()[i] ;
    if (x > old_lower[j]) setColLower(j,x) ; }

  for (i = 0 ; i < n2 ; i++)
  { int j = new_upper.getIndices()[i] ;
    double x = new_upper.getElements()[i] ;
    if (x < old_upper[j]) setColUpper(j,x) ; }
}


void ODSI::deleteCols (int count, const int* cols)

{ if (count <= 0) return ;
/*
  Sort the indices and do the deletions, one by one.
*/
  vector<int> lclcols = vector<int>(&cols[0],&cols[count]) ;

  if (count > 1) std::sort(lclcols.begin(),lclcols.end()) ;

  for (int k = 0 ; k < count ; k++)
  { int j = idx(lclcols[k]) ;
    bool r = consys_delcol(consys, j) ;
    assert(r) ; }
/*
  Now, see if there's an active basis. If so, check that all the variables to
  be deleted are nonbasic. If they are, we can delete them from activeBasis
  and still guarantee a valid basis. If not, throw away activeBasis.
*/
  if (activeBasis)
  { bool allnonbasic = true ;
    OsiDylpWarmStartBasis *odwsb =
      dynamic_cast<OsiDylpWarmStartBasis *>(activeBasis) ;
    for (int k = count-1 ; k >= 0 ; k--)
    { int j = lclcols[k] ;
      if (odwsb->getStructStatus(j) == CoinWarmStartBasis::basic)
      { allnonbasic = false ;
        break ; } }
    if (allnonbasic == true)
    { odwsb->deleteColumns(count,cols) ;
      activeIsModified = true ;
      resolveOptions->forcewarm = true ; }
    else
    { delete activeBasis ;
      activeBasis = 0 ;
      activeIsModified = false ;
      resolveOptions->forcecold = true ; } }

  destruct_cache() ; }


void ODSI::setColSolution (const double* solution)

{ int n = getNumCols() ;
/*
  You can't set what you don't have. Otherwise, replace the current solution.
*/
  if (n == 0) return ;

  if (_col_x) delete[] _col_x ;
  _col_x = new double[n] ;
  
  COPY_VEC(double,solution,_col_x,n) ;
  calc_objval() ;

  return ; }


void ODSI::setRowPrice (const double* price)

{ int m = getNumRows() ;

/*
  You can't set what you don't have. Otherwise, replace the current solution.
*/
  if (m == 0) return ;

  if (_row_price) delete[] _row_price ;
  _row_price = new double[m] ;
  
  COPY_VEC(double,price,_row_price,m) ;

  return ; }




#ifndef _MSC_VER



void ODSI::assert_same (double d1, double d2, bool exact)

{ if (d1 == d2) return ;

  assert(!exact && finite(d1) && finite(d2)) ;

  static const double epsilon = 1.e-10 ;
  double tol = std::max(fabs(d1),fabs(d2))+1 ;
  double diff = fabs(d1 - d2) ;

  assert(diff <= tol*epsilon) ; }


template<class T>
  void ODSI::assert_same (const T& t1, const T& t2, bool exact)


{ if (&t1 == &t2) return ;
  assert(memcmp(&t1,&t2,sizeof(T)) == 0) ; }


template<class T>
  void ODSI::assert_same (const T* t1, const T* t2, int n, bool exact)

{ if (t1 == t2) return ;
  for (int i=0 ; i<n ; i++) assert_same(t1[i],t2[i],exact) ; }


void ODSI::assert_same (const basis_struct& b1, const basis_struct& b2,
                        bool exact)
{ assert(exact) ;                // be explicit
  if (&b1 == &b2) return ;
  assert(b1.len == b2.len) ;

  int size = b1.len*sizeof(basisel_struct) ;
  assert(memcmp(inv_vec<basisel_struct>(b1.el),
                inv_vec<basisel_struct>(b2.el),size) == 0) ; }


void ODSI::assert_same (const conbnd_struct& c1, const conbnd_struct& c2,
                        bool exact)

{ assert(exact) ;                // be explicit
  if (&c1 == &c2) return ;

  //-- consys.h::conbnd_struct
  assert(c1.revs == c2.revs) ;
  assert(c1.inf == c2.inf) ;
  assert(c1.bnd == c2.bnd) ; }


void ODSI::assert_same (const consys_struct& c1, const consys_struct& c2,
                        bool exact)

{ if (&c1 == &c2) return ;
/*
  Simple fields: flags, counts, indices.
*/
  assert(c1.parts == c2.parts) ;
  assert(c1.opts == c2.opts) ;
  assert(c1.varcnt == c2.varcnt) ;
  assert(c1.archvcnt == c2.archvcnt) ;
  assert(c1.logvcnt == c2.logvcnt) ;
  assert(c1.intvcnt == c2.intvcnt) ;
  assert(c1.binvcnt == c2.binvcnt) ;
  assert(c1.maxcollen == c2.maxcollen) ;
  assert(!exact || (c1.maxcolndx == c2.maxcolndx)) ;
  assert(c1.concnt == c2.concnt) ;
  assert(c1.archccnt == c2.archccnt) ;
  assert(c1.cutccnt == c2.cutccnt) ;
  assert(c1.maxrowlen == c2.maxrowlen) ;
  assert(!exact || (c1.maxrowndx == c2.maxrowndx)) ;
  assert(c1.objndx == c2.objndx) ;
  assert(c1.xzndx == c2.xzndx) ;
/*
  Associated vectors.
*/
  int var_count = ODSI::idx(c1.varcnt) ;
  int con_count = ODSI::idx(c1.concnt) ;

  assert_same(c1.obj, c2.obj, var_count, exact) ;
  assert_same(c1.vtyp, c2.vtyp, var_count, true) ;
  assert_same(c1.vub, c2.vub, var_count, exact) ;
  assert_same(c1.vlb, c2.vlb, var_count, exact) ;
  assert_same(c1.colscale, c2.colscale, var_count, exact) ;

  assert_same(c1.rhs, c2.rhs, con_count, exact) ;
  assert_same(c1.rhslow, c2.rhslow, con_count, exact) ;
  assert_same(c1.ctyp, c2.ctyp, con_count, true) ;
  assert_same(c1.cub, c2.cub, con_count, exact) ;
  assert_same(c1.clb, c2.clb, con_count, exact) ;
  assert_same(c1.rowscale, c2.rowscale, con_count, exact) ;
/*
  These can differ and have no mathematical effect. Test them only if the
  client has requested exact equality.
*/
  if (exact)
  { assert(c1.nme == c2.nme) ;
    assert(c1.colsze == c2.colsze) ;
    assert(c1.rowsze == c2.rowsze) ;
    assert(c1.objnme == c2.objnme) ; } }


void ODSI::assert_same (const lpprob_struct& l1, const lpprob_struct& l2,
                        bool exact)

{ if (&l1 == &l2) return ;
/*
  Simple fields: flags, LP phase & return code, objective value, iteration
  count.
*/
  assert(l1.ctlopts == l2.ctlopts) ;
  assert(l1.phase == l2.phase) ;
  assert(l1.lpret == l2.lpret) ;
  assert_same(l1.obj,l2.obj,exact) ;
  assert(l1.iters == l2.iters) ;
/*
  The allocated capacity can differ unless the client requires exact
  equivalence.
*/
  if (exact)
  { assert(l1.colsze == l2.colsze) ;
    assert(l1.rowsze == l2.rowsze) ; }
/*
  Vectors: x, y, status
*/
  int min_col = idx(std::min(l1.colsze, l2.colsze)) ;
  int min_row = idx(std::min(l1.rowsze, l2.rowsze)) ;

  assert_same(l1.status,l2.status,min_col,exact) ;
  assert_same(l1.x,l2.x,min_row, exact) ;
  assert_same(l1.y,l2.y,min_row, exact) ;
/*
  Complex structures: the constraint system and basis.
*/
  assert_same(*l1.consys,*l2.consys,exact) ;
  assert_same(*l1.basis,*l2.basis,exact) ; }


void ODSI::assert_same (const OsiDylpSolverInterface& o1,
                        const OsiDylpSolverInterface& o2, bool exact)

{ 
/*
  The same chunk of memory?
*/
  if (&o1 == &o2) return ;
/*
  These three fields are static. If they don't match, we're really
  confused.
*/
  assert(o1.reference_count == o2.reference_count) ;
  assert(o1.basis_ready == o2.basis_ready) ;
  assert(o1.dylp_owner == o2.dylp_owner) ;
/*
  Test the rest of the simple data fields.
*/
  assert(o1.local_logchn == o2.local_logchn) ;
  assert(o1.initial_gtxecho == o2.initial_gtxecho) ;
  assert(o1.resolve_gtxecho == o2.resolve_gtxecho) ;
  assert(o1.lp_retval == o2.lp_retval) ;
  assert(o1.obj_sense == o2.obj_sense) ;
  assert(o1.solvername == o2.solvername) ;
  assert(o1.mps_debug == o2.mps_debug) ;
/*
  Options, tolerances, and statistics should be byte-for-byte identical.
*/
  assert_same(*o1.initialSolveOptions, *o2.initialSolveOptions, true) ;
  assert_same(*o1.resolveOptions, *o2.resolveOptions, true) ;
  assert_same(*o1.statistics, *o2.statistics, true) ;
  assert_same(*o1.tolerances, *o2.tolerances, true) ;
/*
  The info_ array. The best we can do here is entry by entry equality.
*/
  assert_same(o1.info_,o2.info_,OsiLastHintParam,exact) ;
/*
  Now the more complicated structures. If the lpprob exists, it should
  reference a constraint system, and it should be the same constraint system
  as the consys field. The call to assert_same will check both the lpprob and
  the consys structures.
*/
  if (o1.lpprob)
  { assert(o2.lpprob) ;
    assert(o1.lpprob->consys && (o1.consys == o1.lpprob->consys)) ;
    assert(o2.lpprob->consys && (o2.consys == o2.lpprob->consys)) ;
    assert_same(*o1.lpprob,*o2.lpprob,exact) ; }
  else
  { assert(!o2.lpprob) ;
    assert_same(*o1.consys,*o2.consys,exact) ; }
/*
  A belt & suspenders check: retrieve the matrix as a COINPackedMatrix and
  use the COIN isEquivalent routine to check for equality.
*/
  const CoinPackedMatrix* m1 = o1.getMatrixByCol() ;
  const CoinPackedMatrix* m2 = o2.getMatrixByCol() ;
  if (m1)
  { assert(m2 || m1->isEquivalent(*m2)) ; }
  else
  { assert(!m2) ; }
}


#endif /* !_MSC_VER */



/*
  Problem specification routines: MPS and control files, and loading a problem
  from OSI data structures.
*/


lpopts_struct* main_lpopts ;     // just for cmdint.c::process_cmds
lptols_struct* main_lptols ;     // just for cmdint.c::process_cmds



std::string ODSI::make_filename (const char *filename,
                                 const char *ext1, const char *ext2)

{ std::string basename(filename) ;
  std::string ext1str(ext1), ext2str(ext2) ;

/*
  Extensions must start with ".", so fix if necessary.
*/
  if (ext1 && strlen(ext1) > 0 && ext1[0] != '.')
    ext1str.insert(ext1str.begin(),'.') ;
  if (ext2 && strlen(ext2) > 0 && ext2[0] != '.')
    ext2str.insert(ext2str.begin(),'.') ;
/*
  Strip ext1 from filename, if it exists, leaving the base name.
*/
  if (ext1 && strlen(ext1) > 0)
  { string tmpname(filename) ;
    string::size_type ext1pos = tmpname.rfind(ext1str) ;
    if (ext1pos < tmpname.npos)
    { basename = tmpname.substr(0,ext1pos) ; } }
/*
  Add ext2, if it exists.
*/
  if (ext2 && strlen(ext2) > 0) basename += ext2str ;

  return (basename) ; }








int ODSI::readMps (const char* basename, const char* ext)

/*
  This routine uses the COIN MPS input code, which allows me to sidestep
  the issues that dylp's MPS reader has when confronted with an MPS file
  that requires fixed-field processing or contains a constant objective
  function offset.

  Parameters:
    basename:   base name for mps file
    ext:        file extension; defaults to "mps"

  Returns: -1 if the MPS file could not be opened, otherwise the number of
           errors encountered while reading the file.
*/

{ int errcnt ;
  CoinMpsIO mps ;
  CoinMessageHandler *mps_handler = mps.messageHandler() ;

  if (mps_debug)
  { mps_handler->setLogLevel(handler_->logLevel()) ; }
  else
  { mps_handler->setLogLevel(0) ; }

/*
  See if there's an associated .spc file and read it first, if present.
*/
  std::string filename = make_filename(basename,ext,"spc") ;
  dylp_controlfile(filename.c_str(),true,false) ;
/*
  Make sure the MPS reader has the same idea of infinity as dylp, then
  attempt to read the MPS file.
*/
  mps.setInfinity(DYLP_INFINITY) ;
  filename = make_filename(basename,ext,ext) ;
  errcnt = mps.readMps(filename.c_str(),0) ;
  handler_->message(ODSI_MPSFILEIO,messages_)
    << filename << "read" << errcnt << CoinMessageEol ;
  if (errcnt != 0) return (errcnt) ;
/*
  Load the problem and build a worst-case solution. To take full advantage
  of information in the MPS file, it's easiest to pass the mps object to
  load_problem.
*/
  load_problem(mps) ;
/*
  Done!
*/
  return (0) ; }


void ODSI::writeMps (const char *basename,
                     const char *ext, double objsense) const

{ string filename = make_filename(basename,ext,ext) ;
  
  CoinMpsIO mps ;
  CoinMessageHandler *mps_handler = mps.messageHandler() ;

  if (mps_debug)
  { mps_handler->setLogLevel(handler_->logLevel()) ; }
  else
  { mps_handler->setLogLevel(0) ; }

  double *outputobj ;
  const double *const solverobj = getObjCoefficients() ;
  
  int n = getNumCols(),
      m = getNumRows() ;

  if (objsense == getObjSense())
  { outputobj = const_cast<double *>(solverobj) ; }
  else
  { outputobj = new double[n] ;
    std::transform(solverobj,solverobj+n,outputobj,std::negate<double>()) ; }

  char *vartyp = new char[n] ;
  typedef char *charp ;
  char **colnames = new charp[n],
       **rownames = new charp[m] ;
  int i,j ;

  for (j = 0 ; j < n ; j++) vartyp[j] = isInteger(j) ;

  for (i = 0 ; i < m ; i++)
    rownames[i] = consys_nme(consys,'c',idx(i),false,0) ;
  
  for (j = 0 ; j < n ; j++)
    colnames[j] = consys_nme(consys,'v',idx(j),false,0) ;

  mps.setMpsData(*getMatrixByRow(),DYLP_INFINITY,
                 getColLower(),getColUpper(),outputobj,vartyp,
                 getRowLower(),getRowUpper(),colnames,rownames) ;
/*
  We really need to work on symbolic names for these magic numbers.
*/
  int errcnt = mps.writeMps(filename.c_str(),0,4,2) ;
  handler_->message(ODSI_MPSFILEIO,messages_)
    << filename << "read" << errcnt << CoinMessageEol ;
/*
  Free up the arrays we allocated.
*/
  delete[] vartyp ;
  delete[] colnames ;
  delete[] rownames ;
  if (outputobj != solverobj) delete[] outputobj ;

  return ; }


inline void
ODSI::assignProblem (CoinPackedMatrix*& matrix,
                     double*& col_lower, double*& col_upper, double*& obj, 
                     double*& row_lower, double*& row_upper)

{ loadProblem(*matrix,col_lower,col_upper,obj,row_lower,row_upper) ;
  delete matrix ; matrix = 0 ;
  delete [] col_lower ; col_lower = 0 ;
  delete [] col_upper ; col_upper = 0 ;
  delete [] obj ; obj = 0 ;
  delete [] row_lower ; row_lower = 0 ;
  delete [] row_upper ; row_upper = 0 ;
}


inline void
ODSI::assignProblem (CoinPackedMatrix*& matrix, 
                     double*& lower, double*& upper, double*& obj, 
                     char*& sense, double*& rhs, double*& range)

{ loadProblem(*matrix, lower, upper, obj, sense, rhs, range) ;
  delete matrix ; matrix = 0 ;
  delete [] lower ; lower = 0 ;
  delete [] upper ; upper = 0 ;
  delete [] obj ; obj = 0 ;
  delete [] sense ; sense = 0 ;
  delete [] rhs ; rhs = 0 ;
  delete [] range ; range = 0 ;
}


inline void
ODSI::loadProblem (const CoinPackedMatrix& matrix,
                   const double* col_lower, const double* col_upper,
                   const double* obj,
                   const double* row_lower, const double* row_upper)

{ int rowcnt = matrix.getNumRows() ;
  double *rhs = new double[rowcnt] ;
  double *rhslow = new double[rowcnt] ;
  contyp_enum *ctyp = new contyp_enum[rowcnt] ;

  gen_rowparms(rowcnt,rhs,rhslow,ctyp,row_lower,row_upper) ;
  load_problem(matrix,col_lower,col_upper,obj,ctyp,rhs,rhslow) ;

  delete[] rhs ;
  delete[] rhslow ;
  delete[] ctyp ; }


inline void
ODSI::loadProblem (const CoinPackedMatrix& matrix, 
                   const double* col_lower, const double* col_upper,
                   const double* obj,
                   const char* sense, const double* rhsin, const double* range)

{ int rowcnt = matrix.getNumRows() ;
  double *rhs = new double[rowcnt] ;
  double *rhslow = new double[rowcnt] ;
  contyp_enum *ctyp = new contyp_enum[rowcnt] ;

  gen_rowparms(rowcnt,rhs,rhslow,ctyp,sense,rhsin,range) ;
  load_problem(matrix,col_lower,col_upper,obj,ctyp,rhs,rhslow) ;

  delete[] rhs ;
  delete[] rhslow ;
  delete[] ctyp ; }

inline void
ODSI::loadProblem (const int colcnt, const int rowcnt,
                   const int *start, const int *index, const double *value,
                   const double* col_lower, const double* col_upper,
                   const double* obj,
                   const double* row_lower, const double* row_upper)

{ double *rhs = new double[rowcnt] ;
  double *rhslow = new double[rowcnt] ;
  contyp_enum *ctyp = new contyp_enum[rowcnt] ;

  gen_rowparms(rowcnt,rhs,rhslow,ctyp,row_lower,row_upper) ;
  load_problem(colcnt,rowcnt,start,index,value,
               col_lower,col_upper,obj,ctyp,rhs,rhslow) ;

  delete[] rhs ;
  delete[] rhslow ;
  delete[] ctyp ; }


inline void
ODSI::loadProblem (const int colcnt, const int rowcnt,
                   const int *start, const int *index, const double *value,
                   const double* col_lower, const double* col_upper,
                   const double* obj,
                   const char* sense, const double* rhsin, const double* range)

{ double *rhs = new double[rowcnt] ;
  double *rhslow = new double[rowcnt] ;
  contyp_enum *ctyp = new contyp_enum[rowcnt] ;

  gen_rowparms(rowcnt,rhs,rhslow,ctyp,sense,rhsin,range) ;
  load_problem(colcnt,rowcnt,start,index,value,
               col_lower,col_upper,obj,ctyp,rhs,rhslow) ;

  delete[] rhs ;
  delete[] rhslow ;
  delete[] ctyp ; }







inline double ODSI::getInfinity () const { return DYLP_INFINITY ; }


bool ODSI::setDblParam (OsiDblParam key, double value)
/*
  Aside from simply setting the relevant ODSI field, the work here is making
  sure the dylp tolerance values are consistent.

  OsiDualTolerance and OsiPrimalTolerance are both feasibility tolerances.
  dylp actually maintains feasibility tolerances scaled off the respective
  zero tolerances. dfeas_scale and pfeas_scale allow the scaled tolerances to
  be relaxed (or tightened, for that matter) with respect to the zero
  tolerances. Hence the best approach to consistency is to set *_scale based
  on the ratio of osi and dylp tolerances. This is, at best, imperfect, as it
  doesn't take into account dylp's internal scaling. Nor can it, really, as
  this information isn't guaranteed to be valid.

  Support for Osi[Dual,Primal]ObjectiveLimit is a facade at the moment. The
  values are stored, and used by is[Dual,Primal]ObjectiveLimitReached, but
  dylp knows nothing of them.
*/
{ if (key == OsiLastDblParam) return (false) ;

  bool retval = OSI::setDblParam(key,value) ;

  switch (key)
  { case OsiDualTolerance:
    { tolerances->dfeas_scale = value/tolerances->cost ;
      break ; }
    case OsiPrimalTolerance:
    { tolerances->pfeas_scale = value/tolerances->zero ;
      break ; }
    case OsiObjOffset:
    { break ; }
    case OsiDualObjectiveLimit:
    case OsiPrimalObjectiveLimit:
    { break ; }
    default:
    { retval = false ;
      break ; } }
  
  return (retval) ; }


bool ODSI::getDblParam (OsiDblParam key, double& value) const
/*
  This simply duplicates OSI::getDblParam. I've kept it here until I decide
  what to do with the parameters Osi[Primal,Dual]ObjectiveLimit.
*/
{ if (key == OsiLastDblParam) return (false) ;

  bool retval ;

  switch (key)
  { case OsiDualTolerance:
    case OsiPrimalTolerance:
    case OsiObjOffset:
    { retval = OSI::getDblParam(key,value) ;
      break ; }
    case OsiDualObjectiveLimit:
    case OsiPrimalObjectiveLimit:
    { retval = OSI::getDblParam(key,value) ;
      break ; }
    default:
    { OSI::getDblParam(key,value) ;
      retval = false ;
      break ; } }

  return (retval) ; }


bool ODSI::setIntParam (OsiIntParam key, int value)
/*
  Dylp's iteration limit is normally enforced on a per-phase basis, with an
  overall limit of 3*(per-phase limit). OsiMaxNumIterationHotStart is handled
  in solveFromHotStart.
*/
{ if (key == OsiLastIntParam) return (false) ;

  bool retval = OSI::setIntParam(key,value) ;

  switch (key)
  { case OsiMaxNumIteration:
    { initialSolveOptions->iterlim = value/3 ;
      resolveOptions->iterlim = initialSolveOptions->iterlim ;
      break ; }
    case OsiMaxNumIterationHotStart:
    { break ; }
    default:
    { retval = false ;
      break ; } }
  
  return (retval) ; }


inline bool ODSI::getIntParam (OsiIntParam key, int& value) const
/*
  This simply duplicates OSI::getIntParam. Retained here in anticipation of
  future extensions.
*/
{ bool retval ;

  if (key == OsiLastIntParam) return (false) ;

  switch (key)
  { case OsiMaxNumIteration:
    case OsiMaxNumIterationHotStart:
    { retval = OSI::getIntParam(key,value) ;
      break ; }
    default:
    { retval = false ;
      break ; } }
  
  return (retval) ; }


bool ODSI::setStrParam (OsiStrParam key, const std::string& value)
/*
  Set string paramaters. Note that OsiSolverName is, by definition, a
  constant. You can set it all you like, however :-).

  The problem name is pushed to the constraint system, if it exists.
*/

{ if (key == OsiLastStrParam) return (false) ;

  bool retval = OSI::setStrParam(key,value) ;

  switch (key)
  { case OsiProbName:
    { if (consys)
      { if (consys->nme) strfree(consys->nme) ;
        consys->nme = stralloc(value.c_str()) ; }
      break ; }
    case OsiSolverName:
    { break ; }
    default:
    { retval = false ;
      break ; } }
  
  return (retval) ; }

bool ODSI::getStrParam (OsiStrParam key, std::string& value) const
/*
  This mostly duplicates OSI::getStrParam. It enforces the notion that you
  can't set the solver name, which the default OSI routine does not.
*/
{ if (key == OsiLastStrParam) return (false) ;

  bool retval = true ;

  switch (key)
  { case OsiProbName:
    { retval = OSI::getStrParam(key,value) ;
      break ; }
    case OsiSolverName:
    { value = solvername ;
      break ; }
    default:
    { retval = false ;
      break ; } }

  return (retval) ; }


void ODSI::unimp_hint (bool dylpSense, bool hintSense,
                       OsiHintStrength hintStrength, const char *msgString)

{ if (dylpSense != hintSense)
  { std::string message("Dylp ") ;
    if (dylpSense == true)
    { message += "cannot disable " ; }
    else
    { message += "does not support " ; }
    message += msgString ;
    if (hintStrength == OsiForceDo)
    { handler_->message(ODSI_UNSUPFORCEDO,messages_)
        << message << CoinMessageEol ;
      throw CoinError(message,"setHintParam","OsiDylpSolverInterface") ; }
    else
    { handler_->message(ODSI_IGNORED,messages_)
        << message << CoinMessageEol ; } }
  
  return ; }


bool ODSI::setHintParam (OsiHintParam key, bool sense,
                         OsiHintStrength strength, void *info)
/*
  OSI provides arrays for the sense and strength. ODSI provides the array for
  info.
*/

{ bool retval = false ;
/*
  Set the hint in the OSI structures. Unlike the other set*Param routines,
  setHintParam will return false for key == OsiLastHintParam. Unfortunately,
  it'll also throw for strength = OsiForceDo, without setting a return value.
  We need to catch that throw.
*/
  try
  { retval = OSI::setHintParam(key,sense,strength) ; }
  catch (CoinError)
  { retval = (strength == OsiForceDo) ; }
    
  if (retval == false) return (false) ;
  info_[key] = info ;
/*
  Did the client say `ignore this'? Who are we to argue.
*/
  if (strength == OsiHintIgnore) return (true) ;
/*
  We have a valid hint which would be impolite to simply ignore. Deal with
  it as best we can.
*/
  switch (key)
  {
/*
  Dylp has no presolve capabilities.
*/
    case OsiDoPresolveInInitial:
    case OsiDoPresolveInResolve:
    { unimp_hint(false,sense,strength,"presolve") ;
      retval = true ;
      break ; }
/*
  Dylp can suppress the dual, but cannot suppress the primal. Requesting
  OsiDoDual* with (true, OsiHint*) is ignored; (true,OsiForceDo) will throw
  an exception. On the other hand, suppressing the dual is generally not a
  good idea, so the hint is ignored at (true,OsiHintTry).
*/
    case OsiDoDualInInitial:
    case OsiDoDualInResolve:
    { unimp_hint(false,sense,strength,"exclusive use of dual simplex") ;
      retval = true ;
      lpopts_struct *opts =
        (key == OsiDoDualInInitial)?initialSolveOptions:resolveOptions ;
      if (sense == false)
      { if (strength >= OsiHintDo) opts->usedual = false ; }
      else
      { opts->usedual = true ; }
      break ; }
/*
  No issues here, dylp can go either way.
*/
    case OsiDoScale:
    { if (sense == false)
      { if (strength >= OsiHintTry) initialSolveOptions->scaling = 0 ; }
      else
      { initialSolveOptions->scaling = 2 ; }
      resolveOptions->scaling = initialSolveOptions->scaling ;
      break ; }
/*
  Dylp knows only one crash, which preferably uses slacks, then architecturals,
  and artificials as a last resort.
*/
    case OsiDoCrash:
    { unimp_hint(true,sense,strength,"basis crash") ;
      retval = true ;
      break ; }
/*
  An unadorned hint changes the level by +/- 1, 2, or 3, depending on sense
  and strength (abusing the equivalence of enums and integers).  If info is
  non-zero, it is taken as a pointer to an integer which is the absolute
  value for the log level. In this case, sense is irrelevant. Note that
  CoinMessageHandler recognizes levels 0 -- 4 as generic levels.
  Conveniently, dylp has a similar set of generic levels. The larger powers
  of 2 (8, 16, 32, etc.) left to trigger specific classes of debugging
  printout. These must be set using info. Currently:
    0x08        CoinMpsIO messages
*/
    case OsiDoReducePrint:
    { int verbosity = messageHandler()->logLevel() ;
      mps_debug = false ;
      if (info)
      { verbosity = *reinterpret_cast<int *>(info) ;
        if (verbosity&0x8) mps_debug = true ; }
      else
      { if (sense == true)
        { verbosity -= strength ;
          verbosity = max(verbosity,0) ; }
        else
        { verbosity += strength ;
          verbosity = min(verbosity,4) ; } } 

      dy_setprintopts(0,initialSolveOptions) ;
      dy_setprintopts(0,resolveOptions) ;
      if (verbosity == 0)
      { initial_gtxecho = false ;
        resolve_gtxecho = false ; }
      else
      { initial_gtxecho = true ;
        resolve_gtxecho = true ; }
      messageHandler()->setLogLevel(verbosity) ;
      dy_setprintopts((verbosity&0x7),initialSolveOptions) ;
      dy_setprintopts((verbosity&0x7),resolveOptions) ;
      retval = true ;
      break ; }
/*
  The OSI spec says that unimplemented options (and, by implication, hints)
  should return false. In the case of a hint, however, we can ignore anything
  except OsiForceDo, so usability says we should anticipate new hints and set
  this up so the solver doesn't break. So return true.
*/
    default:
    { unimp_hint(!sense,sense,strength,"unrecognized hint") ;
      retval = true ;
      break ; } }

  return (retval) ; }


inline bool ODSI::getHintParam (OsiHintParam key, bool &sense,
                                OsiHintStrength &strength, void *&info) const
/*
  OSI provides the arrays for the sense and strength. ODSI provides the array
  for info.
*/

{ bool retval = OSI::getHintParam(key,sense,strength) ;

  if (retval == false) return (false) ;

  info = info_[key] ;

  return (true) ; }






void ODSI::initialSolve ()

{ 
/*
  The minimum requirement is a constraint system. Options and tolerances should
  also be present --- they're established when the ODSI object is created.
*/
  assert(consys && initialSolveOptions && tolerances) ;
/*
  Do we have an lpprob structure yet?
*/
  if (!lpprob) construct_lpprob() ;
/*
  Is the basis package initialised?
*/
  if (basis_ready == false)
  { int count = static_cast<int>((1.5*getNumRows()) + 2*getNumCols()) ;
    dy_initbasis(count,initialSolveOptions->factor,0) ;
    basis_ready = true ; }
/*
  Does some other ODSI object own the solver? If so, detach it.
*/
  if (dylp_owner != 0 && dylp_owner != this) dylp_owner->detach_dylp() ;
/*
  Choose options appropriate for the initial solution of the lp and go to it.
*/
  lpopts_struct lclopts = *initialSolveOptions ;
  lptols_struct lcltols = *tolerances ;
  dy_checkdefaults(consys,&lclopts,&lcltols) ;
  if (isactive(local_logchn)) logchn = local_logchn ;
  gtxecho = initial_gtxecho ;
  lpprob->phase = dyINV ;
  lp_retval = dylp_dolp(lpprob,&lclopts,&lcltols,statistics) ;
  if (!(lp_retval == lpOPTIMAL || lp_retval == lpINFEAS ||
        lp_retval == lpUNBOUNDED))
  { throw CoinError("Call to dylp failed.",
                    "initialSolve","OsiDylpSolverInterface") ; }
/*
  Trash the cached solution. This needs to happen regardless of how the lp
  turned out. It needs to be done ahead of getWarmStart, which will try to
  access reduced costs.
*/
    destruct_col_cache() ;
    destruct_row_cache() ;
/*
  Tidy up. If all went well, indicate this object owns the solver, set the
  objective, and set the active basis. If we've failed, do the opposite.
  Note that we have to remove any current active basis first, or getWarmStart
  will hand it back to us.

  Any of lpOPTIMAL, lpINFEAS, or lpUNBOUNDED can (in theory) be hot started
  after allowable problem modifications, hence can be flagged lpctlDYVALID.
  dylp overloads lpprob->obj with the index of the unbounded variable when
  returning lpUNBOUNDED, so we need to fake the objective.
*/
  delete activeBasis ;
  activeBasis = 0 ;
  activeIsModified = false ;
  if (flgon(lpprob->ctlopts,lpctlDYVALID))
  { dylp_owner = this ;
    if (lpprob->lpret == lpUNBOUNDED)
    { _objval = -getObjSense()*getInfinity() ; }
    else
    { _objval = getObjSense()*lpprob->obj ; }
    activeBasis = this->getWarmStart() ; }
  else
  { dylp_owner = 0 ; }

  return ; }






inline int ODSI::getIterationCount () const

{ return ((lpprob)?lpprob->iters:0) ; }


inline bool ODSI::isIterationLimitReached () const

{ return (lp_retval == lpITERLIM) ; }


inline bool ODSI::isProvenOptimal () const

{ return (lp_retval == lpOPTIMAL) ; }


inline bool ODSI::isProvenPrimalInfeasible () const

{ return (lp_retval == lpINFEAS) ; }


inline bool ODSI::isProvenDualInfeasible () const

{ return (lp_retval == lpUNBOUNDED) ; }


inline bool ODSI::isAbandoned () const

{ if (lp_retval == lpACCCHK || lp_retval == lpSINGULAR ||
      lp_retval == lpSTALLED || lp_retval == lpNOSPACE ||
      lp_retval == lpLOSTFEAS || lp_retval == lpPUNT || lp_retval == lpFATAL)
    return (true) ;
  else
    return (false) ; }






inline double ODSI::getObjValue () const
/*
  _objval is kept up-to-date and corrected for obj_sense. All we need to do
  here is deal with the offset.
*/
{ double objoffset ;

  getDblParam(OsiObjOffset,objoffset) ;

  return (_objval-objoffset) ; }
  


const double* ODSI::getColSolution () const
/*
  We could have a cached solution, either supplied by the client or a
  previous solution from dylp. If not, we have to build one.  Dylp returns a
  vector x of basic variables, in basis order, and a status vector. To
  assemble a complete primal solution, it's necessary to construct one vector
  that merges the two sources. The result is cached.

  NOTE the caution above about superbasic (vstatSB) variables. Nonbasic free
  variables (vstatNBFR) occur under the same termination conditions, but are
  legitimately zero. We really should return NaN for superbasics, using the
  function <limits>:std::numeric_limits::signaling_NaN(), but it's just not
  worth the hassle. Sun Workshop C++ (Rogue Wave Software) has it, but Gnu C++
  prior to v3 doesn't even have <limits>, and the local v3 installation's
  <limits> points to another, non-existent header.
*/
{ if (_col_x) return _col_x ;

  if (lpprob)
  { assert(lpprob->status && lpprob->x) ;

    int n = getNumCols() ;
    flags statj ;
    _col_x = new double[n] ;

/*
  Walk the status vector, taking values for basic variables from lpprob.x and
  values for nonbasic variables from the appropriate bounds vector.
*/
    for (int j = 0 ; j < n ; j++)
    { statj = lpprob->status[idx(j)] ;
      if (((int) statj) < 0)
      { int k = -((int) statj) ;
        _col_x[j] = lpprob->x[k] ; }
      else
      { switch (statj)
        { case vstatNBLB:
          case vstatNBFX:
          { _col_x[j] = consys->vlb[idx(j)] ;
            break ; }
          case vstatNBUB:
          { _col_x[j] = consys->vub[idx(j)] ;
            break ; }
          case vstatNBFR:
          { _col_x[j] = 0 ;
            break ; }
          case vstatSB:
          { _col_x[j] = -DYLP_INFINITY ;
            break ; } } } } }

  return (_col_x) ; }


const double* ODSI::getRowPrice () const
/*
  Dylp reports dual variables in basis order, so we need to write a vector
  with the duals dropped into position by row index.
*/

{ if (_row_price) return (_row_price) ;

  if (lpprob)
  { assert(lpprob->basis && lpprob->y) ;

    int m = getNumRows() ;
    _row_price = new double[m] ;
    basis_struct* basis = lpprob->basis ;

    memset(_row_price,0,m*sizeof(double)) ;

    for (int k = 1 ; k <= basis->len ; k++)
    { int i = inv(basis->el[k].cndx) ;
      _row_price[i] = lpprob->y[k]*getObjSense() ; } }

  return (_row_price) ; }


const double *ODSI::getRowActivity () const
/*
  This routine calculates the value of the left-hand-side of a constraint. The
  algorithm is straightforward: Retrieve the primal solution, then walk the
  variables, adding the contribution of each non-zero variable.
*/
{ 
/*
  If we have a cached copy, we're done.
*/
  if (_row_lhs) return (_row_lhs) ;
/*
  In order to go further, we'll need a primal solution and some rows.
*/
  int m = getNumRows() ;
  const double *x = getColSolution() ;

  if (m > 0 && x)
/*
  Create a vector to hold ax and clear it to 0.
*/
  { _row_lhs = new double[consys->concnt] ;
    memset(_row_lhs,0,m*sizeof(double)) ;
/*
  Walk the primal solution vector, adding in the contribution from non-zero
  variables.
*/
    pkvec_struct *aj = pkvec_new(m) ;
    for (int j = 0 ; j < consys->varcnt ; j++)
    { if (x[j] != 0)
      { bool r = consys_getcol_pk(consys,idx(j),&aj) ;
        if (!r)
        { delete[] _row_lhs ;
          if (aj) pkvec_free(aj) ;
          return (0) ; }
        for (int l = 0 ; l < aj->cnt ; l++)
        { int i = inv(aj->coeffs[l].ndx) ;
          _row_lhs[i] += x[j]*aj->coeffs[l].val ; } } }
    if (aj) pkvec_free(aj) ;
/*
  Groom the vector to eliminate tiny values.
*/
    for (int i = 0 ; i < consys->concnt ; i++)
      setcleanzero(_row_lhs[i],tolerances->zero) ; }
/*
  Cache the result and return.
*/

  return (_row_lhs) ; }


const double *ODSI::getReducedCost () const
/*
  Calculate the reduced cost as cbar = c - yA.

  It's tempting to want to dive into dylp for this calculation, but there are
  problems: First, some other ODSI object may own the solver. Second, we
  can't be sure where the duals are coming from --- they could be part of a
  solution from dylp, or they could be a vector supplied by the client.
  Third, both the objective coefficients and the duals are subject to sign
  change for maximisation. In the end, it seems best to rely on
  getObjCoefficients and getRowPrice to provide the objective and duals,
*/

{ 
/*
  Check for a cached value.
*/
  if (_col_cbar) return (_col_cbar) ;
/*
  Do we have columns to to work with? If yes, then we have a valid constraint
  system data structure. Initialize cbar with the objective coefficients.
*/
  int n = getNumCols() ;
  if (n == 0) return (0) ;
  _col_cbar = new double[n] ;
  COPY_VEC(double,getObjCoefficients(),_col_cbar,n) ;  
/*
  Do we have rows? Dual variables? If not, we're done. The presence of rows
  does not necessarily mean we have valid duals.
*/
  int m = getNumRows() ;
  const double *y = getRowPrice() ;
  if (!y) return (_col_cbar) ;
/*
  We have a constraint system, objective coefficients, and dual variables.
  Grind out the calculation, row by row.
*/
  pkvec_struct *ai = pkvec_new(n) ;
  for (int i = 0 ; i < m ; i++)
  { if (y[i] != 0)
    { bool r = consys_getrow_pk(consys,idx(i),&ai) ;
      if (!r)
      { delete[] _col_cbar ;
        _col_cbar = 0 ;
        if (ai) pkvec_free(ai) ;
        return (0) ; }
      for (int l = 0 ; l < ai->cnt ; l++)
      { int j = inv(ai->coeffs[l].ndx) ;
        _col_cbar[j] -= y[i]*ai->coeffs[l].val ; } } }
  if (ai) pkvec_free(ai) ;
/*
  Groom the vector to eliminate tiny values and we're done.
*/
  for (int j = 0 ; j < n ; j++) setcleanzero(_col_cbar[j],tolerances->cost) ;

  return (_col_cbar) ; }






inline bool ODSI::isBinary (int i) const

{ if (!consys || i < 0 || i > consys->varcnt-1)
    return (false) ;
  else
  if (!consys->vtyp)
    return (false) ;
  else
    return (consys->vtyp[idx(i)] == vartypBIN) ; }


inline bool ODSI::isIntegerNonBinary (int i) const

{ if (!consys || i < 0 || i > consys->varcnt-1)
    return (false) ;
  else
  if (!consys->vtyp)
    return (false) ;
  else
    return (consys->vtyp[idx(i)] == vartypINT) ; }


inline bool ODSI::isInteger (int i) const

{ if (!consys || i < 0 || i > consys->varcnt-1)
    return (false) ;
  else
  if (!consys->vtyp)
    return (false) ;
  else
    return (consys->vtyp[idx(i)] == vartypINT ||
            consys->vtyp[idx(i)] == vartypBIN) ; }



inline bool ODSI::isContinuous (int i) const

{ if (!consys || i < 0 || i > consys->varcnt-1)
    return (false) ;
  else
  if (!consys->vtyp)
    return (true) ;
  else
    return (consys->vtyp[idx(i)] == vartypCON) ; }


inline const double* ODSI::getColLower () const

{ if (!consys || !consys->vlb)
    return (0) ;
  else
    return (INV_VEC(double,consys->vlb)) ; }


inline const double* ODSI::getColUpper () const

{ if (!consys || !consys->vub)
    return (0) ;
  else
    return (INV_VEC(double,consys->vub)) ; }


inline int ODSI::getNumCols () const

{ if (!consys)
    return (0) ;
  else
    return (consys->varcnt) ; }


inline int ODSI::getNumElements () const

{ if (!consys)
    return (0) ; 
  else
    return (consys->mtx.coeffcnt) ; }


inline int ODSI::getNumRows () const

{ if (!consys)
    return (0) ;
  else
    return (consys->concnt) ; }


inline const double* ODSI::getObjCoefficients () const

{ if (!consys || !consys->obj) return (0) ;

  if (_col_obj) return _col_obj ;

  int n = getNumCols() ;
  _col_obj = new double[n] ;
  double *obj_0base = INV_VEC(double,consys->obj) ;

  if (getObjSense() < 0)
  { std::transform(obj_0base,obj_0base+n,_col_obj,std::negate<double>()) ; }
  else
  { COPY_VEC(double,obj_0base,_col_obj,n) ; }

  return (_col_obj) ; }


inline double ODSI::getObjSense () const

{ return (obj_sense) ; }


const CoinPackedMatrix* ODSI::getMatrixByRow () const

{ if (!consys) return 0 ;
  if (_matrix_by_row) return _matrix_by_row ;

  _matrix_by_row = new CoinPackedMatrix ;
  _matrix_by_row->reverseOrderedCopyOf(*getMatrixByCol()) ;

  return _matrix_by_row ; }


const CoinPackedMatrix* ODSI::getMatrixByCol () const

{ if (!consys) return 0 ;
  if (_matrix_by_col) return _matrix_by_col ;

/*
  Get the column and coefficient counts and create the OSI core vectors.
*/
  int col_count = getNumCols() ;
  assert(col_count > 0) ;
  int coeff_count = consys->mtx.coeffcnt ;
  assert(coeff_count > 0) ;

  int* start = new int[col_count+1] ;
  int* len = new int[col_count] ;
  double* values = new double[coeff_count] ;
  int* indices = new int[coeff_count] ;
  CoinPackedMatrix* matrix = new CoinPackedMatrix ;
/*
  Scan out the coefficients from the consys_struct column by column.
*/
  colhdr_struct** cols = consys->mtx.cols ;
  int coeff_ndx = 0 ;
  for (int i = 0 ; i < col_count ; i++)
  { start[i] = coeff_ndx ;
    len[i] = cols[idx(i)]->len ;

    coeff_struct* c = cols[idx(i)]->coeffs ;
    for (int j = 0 ; j < len[i] ; j++, c = c->colnxt, coeff_ndx++)
    { values[coeff_ndx] = c->val ;
      indices[coeff_ndx] = inv(c->rowhdr->ndx) ; }
    assert(c == 0) ; }
  assert(coeff_ndx == coeff_count) ;
/*
  Create the OSI sparse matrix structure.
*/
  start[col_count] = coeff_ndx ;
  int row_count = getNumRows() ;
  matrix->assignMatrix(true, row_count, col_count, coeff_count, 
    values, indices, start, len) ;
/*
  Make a cached copy and return a pointer to it.
*/
  _matrix_by_col = matrix ;
  return matrix ; }


const double* ODSI::getRowLower () const

{ if (!consys) return (0) ;
  if (_row_lower) return _row_lower ;

  int n = getNumRows() ;
  double* lower = new double[n] ;

  for (int i = 0 ; i < n ; i++)
  { contyp_enum ctypi = consys->ctyp[idx(i)] ;
    switch (ctypi)
    { case contypEQ:
      case contypGE:
      { lower[i] = consys->rhs[idx(i)] ;
        break ; }
      case contypRNG:
      { lower[i] = consys->rhslow[idx(i)] ;
        break ; }
      case contypLE:
      case contypNB:
      { lower[i] = -DYLP_INFINITY ;
        break ; }
      default:
      { assert(0) ; } } }

  _row_lower = lower ;
  return (lower) ; }


const double* ODSI::getRowUpper () const

{ if (!consys) return (0) ;
  if (_row_upper) return _row_upper ;

  int n = getNumRows() ;
  double* upper = new double[n] ;

  for (int i = 0 ; i < n ; i++)
  { if (consys->ctyp[idx(i)] == contypGE || consys->ctyp[idx(i)] == contypNB)
      upper[i] = DYLP_INFINITY ;
    else
      upper[i] = consys->rhs[idx(i)] ; }

  _row_upper = upper ;
  return upper ; }


const char* ODSI::getRowSense () const

{ if (!consys) return (0) ;
  if (_row_sense) return _row_sense ;
  
  int n = getNumRows() ;
  char* sense = new char[n] ;
  const contyp_enum* ctyp = INV_VEC(contyp_enum,consys->ctyp) ;

  std::transform(ctyp,ctyp+n,sense,type_to_sense) ;

  _row_sense = sense ;
  return sense ; }


const double* ODSI::getRowRange () const

{ if (!consys) return (0) ;
  if (_row_range) return _row_range ;
  
  int n = getNumRows() ;  
  double* range = new double[n] ;
  const double* lower = getRowLower() ;
  const double* upper = getRowUpper() ;
  const char* sense = getRowSense() ;

  for (int i=0 ; i<n ; i++)
    if (sense[i] != 'R')
      range[i] = 0.0 ;
    else
      range[i] = upper[i] - lower[i] ;
  
  _row_range = range ;
  return range ; }


const double* ODSI::getRightHandSide () const

{ if (!consys) return (0) ;
  if (_row_rhs) return _row_rhs ;
  
  int n = getNumRows() ;  
  double* rhs = new double[n] ;
  const double* lower = getRowLower() ;
  const double* upper = getRowUpper() ;
  const char* sense = getRowSense() ;

  for (int i = 0 ; i < n ; i++)
  { switch (sense[i])
    { case 'N':
      { rhs[i] = 0.0 ;
        break ; }
      case 'G':
      { rhs[i] = lower[i] ;
        break ; }
      case 'E':
      { rhs[i] = upper[i] ;
        break ; }
      case 'L':
      { rhs[i] = upper[i] ;
        break ; }
      case 'R':
      { rhs[i] = upper[i] ;
        break ; }
      default:
      { assert(0) ; } } }

  _row_rhs = rhs ;
  return (rhs) ; }







bool ODSI::isDualObjectiveLimitReached () const

{ double objlim ;
  double objval = getObjValue() ;

  getDblParam(OsiDualObjectiveLimit,objlim) ;
  objlim *= getObjSense() ;

  if (getObjSense() > 0)
    return (objval > objlim) ;
  else
    return (objval < objlim) ; }


bool ODSI::isPrimalObjectiveLimitReached () const

{ double objlim ;
  double objval = getObjValue() ;

  getDblParam(OsiDualObjectiveLimit,objlim) ;
  objlim *= getObjSense() ;

  if (getObjSense() > 0)
    return (objval < objlim) ;
  else
    return (objval > objlim) ; }


vector<double*> ODSI::getDualRays (int) const

{ throw CoinError("Unimplemented method.",
                  "getDualRays","OsiDylpSolverInterface") ;

  return (vector<double*>(0)) ; }


vector<double*> ODSI::getPrimalRays (int) const

{ throw CoinError("Unimplemented method.",
                  "getPrimalRays","OsiDylpSolverInterface") ;

  return (vector<double*>(0)) ; }


void ODSI::branchAndBound ()

{ throw CoinError("Unimplemented method.",
                  "branchAndBound","OsiDylpSolverInterface") ;

  return ; }






CoinWarmStart *ODSI::getEmptyWarmStart () const
{ return (dynamic_cast<CoinWarmStart *>(new OsiDylpWarmStartBasis())) ; }

CoinWarmStart* ODSI::getWarmStart () const

/*
  This routine constructs an OsiDylpWarmStartBasis structure from the basis
  returned by dylp.

  Dylp takes the atttitude that in so far as it is possible, logicals should
  be invisible outside the solver. The status of nonbasic logicals is not
  reported, nor does dylp expect to receive it. For completeness, however,
  getWarmStart will synthesize status information for nonbasic logicals.
*/

{ int i,j,k ;
  flags statj ;

/*
  If we have an active basis, return a clone.
*/
  if (activeBasis) { return (activeBasis->clone()) ; }
/*
  Create an empty ODWSB object. If no solution exists, we can return
  immediately.
*/
  OsiDylpWarmStartBasis *wsb = new OsiDylpWarmStartBasis ;
  assert(wsb) ;
  if (!lpprob) return (wsb) ;
/*
  Size the ODWSB object. setSize initialises all status entries to isFree.
  Then grab pointers to the status vectors and the dylp basis vector.
*/
  wsb->setSize(consys->varcnt,consys->concnt) ;

  char *const strucStatus = wsb->getStructuralStatus() ;
  char *const artifStatus = wsb->getArtificialStatus() ;
  char *const conStatus = wsb->getConstraintStatus() ;
  basis_struct *basis = lpprob->basis ;

  if (lpprob->lpret == lpOPTIMAL)
    wsb->setPhase(dyPRIMAL2) ;
  else
    wsb->setPhase(dyPRIMAL1) ;
/*
  Walk the basis and mark the active constraints and basic variables. Basic
  logical variables are entered as the negative of the constraint index.
*/
  for (k = 1 ; k <= basis->len ; k++)
  { i = inv(basis->el[k].cndx) ; 
    setStatus(conStatus,i,CoinWarmStartBasis::atLowerBound) ;
    j = basis->el[k].vndx ;
    if (j < 0)
    { j = inv(-j) ;
      setStatus(artifStatus,j,CoinWarmStartBasis::basic) ; }
    else
    { j = inv(j) ;
      setStatus(strucStatus,j,CoinWarmStartBasis::basic) ; } }
/*
  Next, scan the conStatus array. For each inactive (loose => isFree)
  constraint, indicate that the logical is basic. For active (tight =>
  atLowerBound) constraints, if the corresponding logical isn't already
  marked as basic, use the value of the dual to select one of atLowerBound
  (negative dual) or (for range constraints) atUpperBound (positive dual).
*/
  const double *y = getRowPrice() ;
  for (i = 0 ; i < consys->concnt ; i++)
  { if (getStatus(conStatus,i) == CoinWarmStartBasis::isFree)
    { setStatus(artifStatus,i,CoinWarmStartBasis::basic) ; }
    else
    if (getStatus(artifStatus,i) == CoinWarmStartBasis::isFree)
    { if (y[i] > 0)
      { setStatus(artifStatus,i,CoinWarmStartBasis::atUpperBound) ; }
      else
      { setStatus(artifStatus,i,CoinWarmStartBasis::atLowerBound) ; } } }
/*
  Now scan the status vector and record the status of nonbasic structural
  variables. Some information is lost here --- CoinWarmStartBasis::Status
  doesn't encode NBFX.
*/
  for (j = 1 ; j <= consys->varcnt ; j++)
  { statj = lpprob->status[j] ;
    if (((int) statj) > 0)
    { switch (statj)
      { case vstatNBLB:
        case vstatNBFX:
        { setStatus(strucStatus,inv(j), CoinWarmStartBasis::atLowerBound) ;
          break ; }
        case vstatNBUB:
        { setStatus(strucStatus,inv(j), CoinWarmStartBasis::atUpperBound) ;
          break ; }
        case vstatNBFR:
        { setStatus(strucStatus,inv(j), CoinWarmStartBasis::isFree) ;
          break ; }
        default:
        { delete wsb ;
          wsb = 0 ; } } } }

  return (wsb) ; }


bool ODSI::setWarmStart (const CoinWarmStart *ws)

{ int i,j,k ;
  CoinWarmStartBasis::Status osi_statlogi,osi_statconi,osi_statvarj ;

/*
  Use a dynamic cast to make sure we have an OsiDylpWarmStartBasis. Then check
  for an empty basis.
*/
  const OsiDylpWarmStartBasis *wsb =
                        dynamic_cast<const OsiDylpWarmStartBasis *>(ws) ;
  if (!wsb)
  { handler_->message(ODSI_NOTODWSB,messages_) ;
    return (false) ; }
  int varcnt = wsb->getNumStructural() ;
  int concnt = wsb->getNumArtificial() ;
  if (!(varcnt > 0 && concnt > 0))
  { handler_->message(ODSI_EMPTYODWSB,messages_) ;
    return (false) ; }
/*
  Do we have an lpprob structure yet? If not, construct one.
*/
  assert(consys && consys->vlb && consys->vub) ;
  if (!lpprob) construct_lpprob() ;
  assert(resolveOptions) ;

/*
  Extract the info in the warm start object --- size and status vectors.  The
  number of variables and constraints in the warm start object should match
  the full size of the constraint system. Note that getWarmStart can create
  an empty ODWSB object (see comments with getWarmStart).

  Create a dylp basis_struct and status vector of sufficient size to hold the
  information in the OsiDylpWarmStartBasis. This space may well be freed or
  realloc'd by dylp, so use standard calloc to acquire it.  We'll only use as
  much of the basis as is needed for the active constraints.
*/
  if (!(varcnt == getNumCols() && concnt == getNumRows()))
  { handler_->message(ODSI_ODWSBBADSIZE,messages_) <<
      concnt << varcnt << getNumRows() << getNumCols() ;
    return (false) ; }

  const char *const strucStatus = wsb->getStructuralStatus() ;
  const char *const artifStatus = wsb->getArtificialStatus() ;
  const char *const conStatus = wsb->getConstraintStatus() ;

  flags *status = new flags[idx(varcnt)] ;
  basis_struct basis ;
  basis.el = new basisel_struct[idx(concnt)] ;
/*
  We never use these entries, but we need to initialise them to squash a
  `read from uninitialised memory' error during block copies.
*/
  basis.el[0].cndx = 0 ;
  basis.el[0].vndx = 0 ;
  status[0] = 0 ;
/*
  Walk the constraint status vector and build the set of active constraints.
*/
  int actcons = 0 ;
  for (i = 1 ; i <= concnt ; i++)
  { osi_statconi = getStatus(conStatus,inv(i)) ;
    if (osi_statconi == CoinWarmStartBasis::atLowerBound)
    { actcons++ ;
      basis.el[actcons].cndx = i ;
      basis.el[actcons].vndx = 0 ; } }
  basis.len = actcons ;
/*
  Now walk the structural status vector. For each basic variable, drop it into
  a basis entry and note the position in the status vector. For each nonbasic
  variable, set the proper status flag. We need to check bounds to see if the
  variable should be fixed.
*/
  k = 0 ;
  for (j = 1 ; j <= varcnt ; j++)
  { osi_statvarj = getStatus(strucStatus,inv(j)) ;
    switch (osi_statvarj)
    { case CoinWarmStartBasis::basic:
      { k++ ;
        assert(k <= actcons) ;
        basis.el[k].vndx = j ;
        status[j] = (flags) (-k) ;
        break ; }
      case CoinWarmStartBasis::atLowerBound:
      { if (consys->vlb[j] == consys->vub[j])
        { status[j] = vstatNBFX ; }
        else
        { status[j] = vstatNBLB ; }
        break ; }
      case CoinWarmStartBasis::atUpperBound:
      { status[j] = vstatNBUB ;
        break ; }
      case CoinWarmStartBasis::isFree:
      { status[j] = vstatNBFR ;
        break ; } } }
/*
  Now we need to finish out the basis by adding the basic logicals.  The
  convention is to represent this with the negative of the index of the
  constraint which spawns the logical. Note that we do this only if the
  constraint is active. When running in fullsys mode, this is the common
  case. If we're running dynamic simplex, then all we're picking up here
  is the occasional logical that's basic at bound.
*/
  for (i = 1 ; i <= concnt ; i++)
  { osi_statlogi = getStatus(artifStatus,inv(i)) ;
    osi_statconi = getStatus(conStatus,inv(i)) ;
    if (osi_statlogi == CoinWarmStartBasis::basic &&
        osi_statconi == CoinWarmStartBasis::atLowerBound)
    { k++ ;
      assert(k <= actcons) ;
      basis.el[k].vndx = -i ; } }
/*
  Now install the new basis in the lpprob_struct. If the present allocated
  capacity is sufficient, we can just copy over the existing information. If
  the capacity is insufficient, we'll free the existing space and allocate
  new. There's also the possibility that this WarmStartBasis came from some
  other object, and there are no vectors here at all. Because dylp relies on
  lpprob->colsze and lpprob->rowsze for the allocated capacity of actvars, x,
  and y, we need to reallocate them too.
*/
  if (lpprob->colsze < varcnt)
  { if (lpprob->status)
    { FREE(lpprob->status) ;
      lpprob->status = 0 ; }
    if (lpprob->actvars)
    { lpprob->actvars =
        (bool *) REALLOC(lpprob->actvars,idx(varcnt)*sizeof(bool)) ; }
    lpprob->colsze = varcnt ; }
  if (!lpprob->status)
  { lpprob->status = (flags *) CALLOC(idx(lpprob->colsze),sizeof(flags)) ; }
  if (!lpprob->actvars)
  { lpprob->actvars = (bool *) CALLOC(idx(lpprob->colsze),sizeof(bool)) ; }

  if (lpprob->rowsze < actcons)
  { if (lpprob->x)
    { lpprob->x =
        (double *) REALLOC(lpprob->x,idx(actcons)*sizeof(double)) ; }
    if (lpprob->y)
    { lpprob->y =
        (double *) REALLOC(lpprob->y,idx(actcons)*sizeof(double)) ; }
    if (lpprob->basis && lpprob->basis->el)
    { FREE(lpprob->basis->el) ;
      lpprob->basis->el = 0 ; }
    lpprob->rowsze = actcons ; }
  if (!lpprob->x)
  { lpprob->x = (double *) CALLOC(idx(lpprob->rowsze),sizeof(double)) ; }
  if (!lpprob->y)
  { lpprob->y = (double *) CALLOC(idx(lpprob->rowsze),sizeof(double)) ; }
  if (!lpprob->basis)
  { lpprob->basis =  (basis_struct *) CALLOC(1,sizeof(basis_struct)) ; }
  if (!lpprob->basis->el)
  { lpprob->basis->el =
      (basisel_struct *) CALLOC(idx(lpprob->rowsze),sizeof(basisel_struct)) ; }
/*
  Whew. The actual copy is dead easy.
*/
  copy_basis(&basis,lpprob->basis) ;
  COPY_VEC(flags,status,lpprob->status,idx(varcnt)) ;
  delete[] basis.el ;
  delete[] status ;
/*
  And we're done. The new status and basis are installed in lpprob, and the
  x, y, and actvars arrays have been resized if needed. Set the phase, signal
  a warm start, and then we're done.
*/
  lpprob->phase = wsb->getPhase() ;
  resolveOptions->forcecold = false ;
  resolveOptions->forcewarm = true ;
/*
  One last thing --- make a copy of wsb to be the active basis (that is, unless
  wsb is the active basis, and we're simply installing it in lpprob).
*/
  if (wsb != static_cast<const CoinWarmStart *>(activeBasis) )
  { delete activeBasis ;
    activeBasis = wsb->clone() ;
    activeIsModified = false ; }
  
  return (true) ; }


void ODSI::resolve ()
/*
  Grossly oversimplified, this routine calls dylp after making sure that the
  forcecold option is turned off.  But we do need to make sure the solver is
  ours to use, and there's some state to maintain.

  If we're reoptimising, then the basis should be ready and we should have
  warm start information. For our purposes here, that boils down to the
  presence of an activeBasis object (created by setWarmStart() or a recent call
  to the solver).

  Note that we don't actually force a warm start unless we install the active
  basis here. Similarly, setWarmStart will force a warm start. But if we have
  an unmodified active basis, dylp should be able to manage a hot start.
*/
{ assert(lpprob && lpprob->basis && lpprob->status && consys &&
         resolveOptions && tolerances && statistics) ;

  if (dylp_owner != 0 && dylp_owner != this) dylp_owner->detach_dylp() ;
/*
  Is the basis package initialised?
*/
  if (basis_ready == false)
  { int count = static_cast<int>((1.5*getNumRows()) + 2*getNumCols()) ;
    dy_initbasis(count,initialSolveOptions->factor,0) ;
    basis_ready = true ; }
/*
  We can hope that activeBasis is already installed in the lpprob, but there
  are several circumstances to trap for here:
    * The user has been playing with cuts since the last call to dylp, but
      kept within the rules for modifying the active basis. In this case
      activeIsModified will be true, and we need to install it.
    * The solver has been used by some other ODSI object since the last
      call by this object. In this case, we've just done a detach_dylp()
      and dylp_owner will be null. Again, we need to install activeBasis.
  If we have an active basis, install it. If that fails, remove activeBasis.
*/
  if (activeBasis && (activeIsModified == true || dylp_owner == 0))
  { if (setWarmStart(activeBasis) == false)
    { throw CoinError("Warm start failed --- invalid active basis.",
                      "resolve","OsiDylpSolverInterface") ; }
    delete activeBasis ;
    activeBasis = 0 ;
    activeIsModified = false ;
    resolveOptions->forcewarm = true ; }
/*
  Choose options appropriate for reoptimising and go to it. Phase is not
  crucial here --- dylp will sort it out as it starts.
*/
  lpopts_struct lclopts = *resolveOptions ;
  lptols_struct lcltols = *tolerances ;
  dy_checkdefaults(consys,&lclopts,&lcltols) ;
  assert(lclopts.forcecold == false) ;
  if (isactive(local_logchn)) logchn = local_logchn ;
  gtxecho = resolve_gtxecho ;
  dyphase_enum phase = lpprob->phase ;
  if (!(phase == dyPRIMAL1 || phase == dyPRIMAL2 || phase == dyDUAL))
    lpprob->phase = dyINV ;
  lp_retval = dylp_dolp(lpprob,&lclopts,&lcltols,statistics) ;
  if (!(lp_retval == lpOPTIMAL || lp_retval == lpINFEAS ||
        lp_retval == lpUNBOUNDED))
  { throw CoinError("Call to dylp failed.",
                    "resolve","OsiDylpSolverInterface") ; }
/*
  Trash the cached solution. This needs to happen regardless of how the lp
  turned out. It needs to be done ahead of getWarmStart, which will try to
  access reduced costs.
*/
  destruct_col_cache() ;
  destruct_row_cache() ;
/*
  Tidy up. If all went well, indicate this object owns the solver, set the
  objective, and set the active basis. If we've failed, do the opposite.
  Note that we have to remove any current active basis first, or getWarmStart
  will hand it back to us.

  Any of lpOPTIMAL, lpINFEAS, or lpUNBOUNDED can (in theory) be hot started
  after allowable problem modifications, hence can be flagged lpctlDYVALID.
  dylp overloads lpprob->obj with the index of the unbounded variable when
  returning lpUNBOUNDED, so we need to fake the objective.
*/
  delete activeBasis ;
  activeBasis = 0 ;
  activeIsModified = false ;
  if (flgon(lpprob->ctlopts,lpctlDYVALID))
  { dylp_owner = this ;
    if (lpprob->lpret == lpUNBOUNDED)
    { _objval = -getObjSense()*getInfinity() ; }
    else
    { _objval = getObjSense()*lpprob->obj ; }
    activeBasis = this->getWarmStart() ;
    resolveOptions->forcewarm = false ; }
  else
  { dylp_owner = 0 ; }
  
  return ; }





inline void ODSI::markHotStart ()
{ 
/*
  If we're attempting this, then it should be true that this OSI object
  owns the solver, has successfully solved an lp, and has left dylp's data
  structures intact.
*/
  assert(dylp_owner == this) ;
  assert(lpprob && resolveOptions) ;
  assert(flgon(lpprob->ctlopts,lpctlDYVALID)) ;
/*
  Turn off the forcecold and forcewarm options, allowing a hot start.
*/
  resolveOptions->forcecold = false ;
  resolveOptions->forcewarm = false ;
/*
  Capture a warm start object, in case some other ODSI object uses the solver
  between hot starts by this ODSI object.
*/
  if (hotstart_fallback) delete hotstart_fallback ;
  hotstart_fallback = getWarmStart() ;

  return ; }

void ODSI::solveFromHotStart ()

{ 
/*
  If some other ODSI object has used the solver, then fall back to a warm
  start. Throw an error if there's no warm start object to fall back on.
  (The throw message lies a little, for the benefit of users who haven't
  checked the details of the code.)
*/

  if (dylp_owner != this)
  { if (hotstart_fallback && setWarmStart(hotstart_fallback))
    { resolve() ; }
    else
    { throw CoinError("Hot start failed --- invalid/missing hot start object.",
                      "solveFromHotStart","OsiDylpSolverInterface") ; } }
/*
  If no other ODSI object has used the solver, all we need to do is
  make sure forcecold and forcewarm are off. The basis should be ready.
  Note that dylp does need to know what's changed: any of bounds, rhs &
  rhslow, or objective. The various routines that make these changes
  set the flags.
*/
  else
  { int tmp_iterlim = -1 ;
    int hotlim ;

    assert(lpprob && lpprob->basis && lpprob->status && basis_ready &&
           consys && resolveOptions && tolerances && statistics) ;

    if (isactive(local_logchn)) logchn = local_logchn ;
    gtxecho = resolve_gtxecho ;
/*
  Setting the phase to dyPRIMAL1 is overly cautious, but the only practical
  solution (dylp expects the client to set this, but exposing this would
  violates the generic OSI interface). In any event, dylp will sort it out.
*/
    lpopts_struct lclopts = *resolveOptions ;
    lptols_struct lcltols = *tolerances ;
    dy_checkdefaults(consys,&lclopts,&lcltols) ;
    lpprob->phase = dyPRIMAL1 ;
    assert(lclopts.forcecold == false) ;
    lclopts.forcewarm = false ;
    getIntParam(OsiMaxNumIterationHotStart,hotlim) ;
    if (hotlim > 0)
    { tmp_iterlim = resolveOptions->iterlim ;
      resolveOptions->iterlim = (hotlim/3 > 0)?hotlim/3:1 ; }

    lp_retval = dylp_dolp(lpprob,&lclopts,&lcltols,statistics) ;
    if (!(lp_retval == lpOPTIMAL || lp_retval == lpINFEAS ||
          lp_retval == lpUNBOUNDED))
    { throw CoinError("Call to dylp failed.",
                      "solveFromHotStart","OsiDylpSolverInterface") ; }
/*
  Trash the cached solution. This needs to happen regardless of how the lp
  turned out. It needs to be done ahead of getWarmStart, which will try to
  access reduced costs.
*/
    destruct_col_cache() ;
    destruct_row_cache() ;
/*
  Tidy up. If all went well, indicate this object owns the solver, set the
  objective, and set the active basis. If we've failed, do the opposite.
  Note that we have to remove any current active basis first, or getWarmStart
  will hand it back to us.

  Any of lpOPTIMAL, lpINFEAS, or lpUNBOUNDED can (in theory) be hot started
  after allowable problem modifications, hence can be flagged lpctlDYVALID.
  dylp overloads lpprob->obj with the index of the unbounded variable when
  returning lpUNBOUNDED, so we need to fake the objective.
*/
    delete activeBasis ;
    activeBasis = 0 ;
    activeIsModified = false ;
    if (flgon(lpprob->ctlopts,lpctlDYVALID))
    { if (lpprob->lpret == lpUNBOUNDED)
      { _objval = -getObjSense()*getInfinity() ; }
      else
      { _objval = getObjSense()*lpprob->obj ; }
      activeBasis = this->getWarmStart() ; }
    else
    { dylp_owner = 0 ; }
    if (tmp_iterlim > 0) resolveOptions->iterlim = tmp_iterlim ; }

  return ; }


inline void ODSI::unmarkHotStart ()
/*
  Nothing to do here but to destroy the fall back warm start object.
*/

{ if (hotstart_fallback)
  { delete hotstart_fallback ;
    hotstart_fallback = 0 ; }

  return ; }





void ODSI::activateRowCutDebugger (const char * modelName)

{ delete rowCutDebugger_ ;

  if (dylp_owner && dylp_owner->lpprob &&
      flgon(dylp_owner->lpprob->ctlopts,lpctlDYVALID))
  { CoinWarmStart *ws = dylp_owner->getWarmStart() ;
    ODSI *prev_owner = dylp_owner ;
    prev_owner->detach_dylp() ;
    rowCutDebugger_ = new OsiRowCutDebugger(*this,modelName) ;
    prev_owner->setWarmStart(ws) ;
    prev_owner->resolve() ;
    delete ws ; }
  else
  { rowCutDebugger_ = new OsiRowCutDebugger(*this,modelName) ; }

  return ; }





void ODSI::dylp_controlfile (const char *name,
                             const bool silent, const bool mustexist)

{ if (name == 0 || *name == 0) return ;
  string mode = (mustexist)?"r":"q" ;
  cmdchn = openfile(const_cast<char *>(name),const_cast<char *>(mode.c_str())) ;
  if (!(cmdchn == IOID_INV || cmdchn == IOID_NOSTRM))
  { setmode (cmdchn, 'l') ;  
    main_lpopts = initialSolveOptions ;
    main_lptols = tolerances ;
    bool r = (process_cmds(silent) != 0) ;
    (void) closefile(cmdchn) ;
    assert(r == cmdOK) ;
/*
  Copy user settings into resolveOptions, except for forcecold and fullsys.
*/
    lpopts_struct saveOptions = *resolveOptions ;
    memcpy(resolveOptions,initialSolveOptions,sizeof(lpopts_struct));
    resolveOptions->forcecold = saveOptions.forcecold ;
    resolveOptions->fullsys = saveOptions.fullsys ; }

  cmdchn = IOID_NOSTRM ;

  return ; }
  

void ODSI::dylp_logfile (const char *name, bool echo)

{ if (name == 0 || *name == 0) return ;

  string log = make_filename(name,".mps",".log") ;
  local_logchn = openfile(const_cast<char *>(log.c_str()),
                          const_cast<char *>("w")) ;
  if (local_logchn == IOID_INV) local_logchn = IOID_NOSTRM ;
  initial_gtxecho = echo ;
  resolve_gtxecho = echo ; }


#endif /* COIN_USE_DYLP */

---