/* $Id: CoinPresolvePsdebug.cpp 1560 2012-11-24 00:29:01Z lou $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#include <new>
#include <stdio.h>
#include <math.h>

#include "CoinPresolveMatrix.hpp"
#include "CoinHelperFunctions.hpp"

/*
  \file

  This file contains a number of routines that are useful when doing
  serious debugging but unneeded otherwise. It also contains the methods
  that implement the CoinPresolveMonitor class.

  Presumably if you're deep enough into presolve to need these, you're
  willing to scan the file to see what can be done. See also the Presolve
  Debug Functions module in the doxygen doc'n and CoinPresolvePsdebug.hpp.

  The general approach for the matrix consistency routines is that the
  routines return void and abort when they find a problem. The routines
  that check the basis and solution complain loudly but do not abort.

  NOTE: The definitions for PRESOLVE_CONSISTENCY and PRESOLVE_DEBUG MUST
	BE CONSISTENT across all CoinPresolve source files AND OsiPresolve
	AND ClpPresolve AND OsiDylpPresolve (assuming any of the latter
	are relevant to your code).  Otherwise, at best you'll get garbage
	output. More likely, you'll get a core dump. Resist the temptation
	to define these constants in individual files. In particular,
	cdone and rdone will NOT be consistently maintained during postsolve.

  Hack away as your needs dictate.
*/


/*
  Integrity checking routines for the (loosely) packed matrices of a
  CoinPresolveMatrix object. Some routines work on column-major and row-major
  reps separately, others do cross-checking.
*/

namespace { // begin unnamed file-local namespace

#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
/*
  Check for duplicate entries in a major vector by walking the vector. For
  each coefficient, use presolve_find_minor1 to search the remainder of the
  major vector for an entry with the same minor index. We don't want to find
  anything.
*/
  
void no_majvec_dups (const char *majdones, const CoinBigIndex *majstrts,
		     const int *minndxs, const int *majlens, int nmaj)

{ for (int maj = 0 ; maj < nmaj ; maj++) 
  { if ((!majdones || majdones[maj]) && majlens[maj] > 0)
    { CoinBigIndex ks = majstrts[maj] ;
      CoinBigIndex ke = ks+majlens[maj] ;
      for (CoinBigIndex k = ks ; k < ke ; k++)
      { 
/*
  Assert we fell off the end of the major vector without finding the entry. 
*/
	PRESOLVEASSERT(presolve_find_minor1(minndxs[k],k+1,
					    ke,minndxs) == ke) ; } } }
  return ; }

/*
  As the name implies: scan for explicit zeros.
*/
void check_majvec_nozeros (const CoinBigIndex *majstrts, const double *majels,
			   const int *majlens, int nmaj)

{ for (int maj = 0 ; maj < nmaj ; maj++) 
  { if (majlens[maj] > 0)
    { CoinBigIndex ks = majstrts[maj] ;
      CoinBigIndex ke = ks+majlens[maj] ;
      for (CoinBigIndex k = ks ; k < ke ; k++)
      { PRESOLVEASSERT(fabs(majels[k]) > ZTOLDP) ; } } }

  return ; }

/*
  Integrity checks for the linked lists that indicate major vector ordering
  in the bulk storage area (minor index and coefficient arrays).
 */
void links_ok (presolvehlink *majlink, int *majstrts, int *majlens, int nmaj)

{ int maj ;

/*
  Confirm link integrity. Vectors of length 0 should not be part of the chain.
*/
  for (maj = 0 ; maj < nmaj ; maj++)
  { int pre = majlink[maj].pre ;
    int suc = majlink[maj].suc ;

    if (majlens[maj] == 0)
    { PRESOLVEASSERT(pre == NO_LINK && suc == NO_LINK) ; }
    if (pre != NO_LINK)
    { PRESOLVEASSERT(0 <= pre && pre <= nmaj) ;
      PRESOLVEASSERT(majlink[pre].suc == maj) ; }
    if (suc != NO_LINK)
    { PRESOLVEASSERT(0 <= suc && suc <= nmaj) ;
      PRESOLVEASSERT(majlink[suc].pre == maj) ; } }
/*
  There must be a first vector.
*/
  for (maj = 0 ; maj < nmaj ; maj++) 
  { if (majlink[maj].pre == NO_LINK)
      break ; }
  PRESOLVEASSERT(nmaj == 0 || maj < nmaj) ;
/*
  The order of the linked list should match the ordering indicated by the
  major vector start & length arrays.
*/
  while (maj != NO_LINK)
  { if (majlink[maj].suc != NO_LINK) 
    { PRESOLVEASSERT(majstrts[maj]+majlens[maj] <=
				majstrts[majlink[maj].suc]) ; }
    maj = majlink[maj].suc ; }

  return ; }


/*
 matrix_consistent checks that an entry is in the column-major representation
 if it is in the row-major representation.  If testvals is non-zero, it also
 checks that their values are the same.

 By doing the appropriate swaps of column- and row-major data structures in
 the parameter list, we can check that an entry is in the row-major
 representation if it's in the column-major representation.

 I can't see any nice way to rename the parameters (majmajstrt? minmajstrt?).

 Original comment:  ``Note that there may be entries in a row that correspond
 to empty columns and vice-versa.'' To which a previous browser had commented
 ``HUH???''. And I agree. -- lh, 040907 --
*/

void matrix_consistent (const CoinBigIndex *mrstrt, const int *hinrow,
			const int *hcol, const double *rowels,
			const CoinBigIndex *mcstrt, const int *hincol,
			const int *hrow, const double *colels,
			int nrows, int testvals,
			const char *ROW, const char *COL)
{
  for (int irow = 0 ; irow < nrows ; irow++) {
    if (hinrow[irow] > 0) {
      const CoinBigIndex krs = mrstrt[irow] ;
      const CoinBigIndex kre = krs+hinrow[irow] ;

      for (CoinBigIndex k = krs ; k < kre ; k++) {
	int jcol = hcol[k] ;
	const CoinBigIndex kcs = mcstrt[jcol] ;
	const CoinBigIndex kce = kcs+hincol[jcol] ;

	CoinBigIndex kk = presolve_find_row1(irow,kcs,kce,hrow) ;
	if (kk == kce) {
	  std::cout
	    << "MATRIX INCONSISTENT:  can't find " << ROW << " " << irow
	    << " in " << COL << " " << jcol << std::endl ;
	  fflush(stdout) ;
	  abort() ;
	}
	if (testvals && colels[kk] != rowels[k]) {
	  std::cout
	    << "MATRIX INCONSISTENT: values differ for " << ROW << " " << irow
	    << " and " << COL << " " << jcol << std::endl ;
	  fflush(stdout) ;
	  abort() ;
	}
      }
    }
  }
}
#endif
} // end unnamed file-local namespace

/*
  Utilizes matrix_consistent to check for equivalence of the column- and
  row-major representations. Checks for presence of coefficients in the
  column-major matrix, given presence in the row-major matrix, then checks
  for presence in the row-major matrix given presence in the column-major
  matrix. If testvals == true (default), the check also tests that the
  coefficients have equal value.
  
  See further comments with matrix_consistent.
*/

# if PRESOLVE_CONSISTENCY
void presolve_consistent(const CoinPresolveMatrix *preObj, bool testvals)
{
  matrix_consistent(preObj->mrstrt_,preObj->hinrow_,preObj->hcol_,
		    preObj->rowels_,
		    preObj->mcstrt_,preObj->hincol_,preObj->hrow_,
		    preObj->colels_,
		    preObj->nrows_,testvals,"row","col") ;
  matrix_consistent(preObj->mcstrt_,preObj->hincol_,preObj->hrow_,
		    preObj->colels_,
		    preObj->mrstrt_,preObj->hinrow_,preObj->hcol_,
		    preObj->rowels_, 
		    preObj->ncols_,testvals,"col","row") ;
}

/*
  Check the column- and/or row-major matrices for duplicates. By default, both
  will be checked.
*/
  
void presolve_no_dups (const CoinPresolveMatrix *preObj,
		       bool doCol, bool doRow)

{
  if (doCol)
  { no_majvec_dups(0,preObj->mcstrt_,preObj->hrow_,
		   preObj->hincol_,preObj->ncols_) ; }
  if (doRow)
  { no_majvec_dups(0,preObj->mrstrt_,preObj->hcol_,
		   preObj->hinrow_,preObj->nrows_) ; }

  return ; }


/*
  As the name implies: scan for explicit zeros. By default, both matrices are
  scanned.
*/
void presolve_no_zeros (const CoinPresolveMatrix *preObj,
			bool doCol, bool doRow)
{
  if (doCol)
  { check_majvec_nozeros(preObj->mcstrt_,preObj->colels_,preObj->hincol_,
			 preObj->ncols_) ; }
  if (doRow)
  { check_majvec_nozeros(preObj->mrstrt_,preObj->rowels_,preObj->hinrow_,
			 preObj->nrows_) ; }

  return ; }

/*
  Lazy check on column lengths. Scan the row index array for the column.
  If the relevant row length in the row-major rep is non-zero, assume we're ok.

  Not advertised in CoinPresolvePsdebug.hpp.
*/
void presolve_hincol_ok(const int *mcstrt, const int *hincol,
	       const int *hinrow,
	       const int *hrow, int ncols)
{
  int jcol;

  for (jcol=0; jcol<ncols; jcol++) 
    if (hincol[jcol] > 0) {
      int kcs = mcstrt[jcol];
      int kce = kcs + hincol[jcol];
      int n=0;
      
      int k;
      for (k=kcs; k<kce; k++) {
	int row = hrow[k];
	if (hinrow[row] > 0)
	  n++;
      }
      if (n != hincol[jcol])
	abort();
    }
}

/*
  Integrity checks for the linked lists that indicate major vector ordering
  in the bulk storage area (minor index and coefficient arrays).
 */
void presolve_links_ok (const CoinPresolveMatrix *preObj,
			bool doCol, bool doRow)
{
  if (doCol)
  { links_ok(preObj->clink_,preObj->mcstrt_,
	     preObj->hincol_,preObj->ncols_) ; }
  if (doRow)
  { links_ok(preObj->rlink_,preObj->mrstrt_,
	     preObj->hinrow_,preObj->nrows_) ; }

  return ; }



/*
  Routines to check a threaded matrix from a CoinPostsolve object.
*/

/*
  Check that the column length agrees with the column thread. There must be
  the correct number of coefficients, and the thread must end with the NO_LINK
  marker.
*/
void presolve_check_threads (const CoinPostsolveMatrix *obj)

{ 

  CoinBigIndex *mcstrt = obj->mcstrt_ ;
  int *hincol = obj->hincol_ ;
  CoinBigIndex *link = obj->link_ ;
  char *cdone = obj->cdone_ ;

  int n = obj->ncols0_ ;

/*
  Scan the columns, checking only the ones that have been processed into the
  constraint matrix.
*/
  for (int j = 0 ; j < n ; j++)
  { if (!cdone[j]) continue ;

    int lenj = hincol[j] ;
    int k ;
    for (k = mcstrt[j] ; k != NO_LINK && lenj > 0 ; k = link[k])
    { assert(k >= 0 && k < obj->maxlink_) ;
      lenj-- ; }

    assert(k == NO_LINK && lenj == 0) ; }

  return ; }

/*
  Check the free list. We're looking for gross corruption here. The notion is
  that the free list plus elements in the matrix should add up to the capacity
  of the bulk store.
*/

void presolve_check_free_list (const CoinPostsolveMatrix *obj, bool chkElemCnt)

{ 

  CoinBigIndex k = obj->free_list_ ;
  CoinBigIndex freeCnt = 0 ;
  CoinBigIndex maxlink = obj->maxlink_ ;
  CoinBigIndex *link = obj->link_ ;
/*
  Redundancy in the data structure. These should always be equal.
*/
  assert(maxlink == obj->bulk0_) ;
/*
  Walk the free list portion of link. We should never point outside the bulk
  store. If we ever come across an entry that's less than 0, it had better be
  NO_LINK, the end marker.
*/
  while (k >= 0)
  { assert(k < maxlink) ;
    freeCnt++ ;
    k = link[k] ; }
  assert(k == NO_LINK) ;
/*
  And a final test: elements in the matrix plus free space should equal the
  size of the bulk area. A good thought, but less than practical. Currently
  postsolve doesn't track the number of elements in the matrix. But you might
  find it useful if you're checking a newly constructed postsolve matrix. Even
  then, you need to make sure nelems_ is correct. In the normal scheme of
  things, this requires that somewhere there's a count of elements. Right now,
  drop_empty_cols_action::presolve does this count, and you can get an accurate
  value from the presolve object. assignPresolveToPostsolve will transfer this
  value. Otherwise you're on your own --- your constructor must somehow find
  this count. Using a standard CoinPackedMatrix is another way to get a count.
*/
  if (chkElemCnt)
  { assert(obj->nelems_+freeCnt == maxlink) ; }


  return ; }
#endif


/*
  Routines to check solution and basis composition.
*/

/*
  CoinPostsolveMatrix

  This routine performs two checks on reduced costs held in rcosts_:
    * The value held in rcosts_ is checked against the status of the
      variable. Errors reported as "Bad rcost"
    * The reduced cost is calculated from scratch and compared to the
      value held in rcosts_. Errors reported as "Inacc rcost"

  Remember that postsolve has a schizophrenic attitude about maximisation. All
  transforms assume minimisation, and that's reflected in the reduced costs we
  see here. And you must load duals and reduced costs with the correct sign for
  minimisation. But, as a small courtesy (and a big inconsistency), postsolve
  will negate objective coefficients for you. Hence the rather odd use of
  maxmin.

  The routine is specific to CoinPostsolveMatrix because the reduced cost
  calculation requires traversal of (threaded) matrix columns.

  NOTE: This routine holds static variables. It will detect when the problem
	size changes and reinitialise. If you use presolve debugging over
	multiple problems and you want to be dead sure of reinitialisation,
	use the call presolve_check_reduced_costs(0), which will reinitialise
	and return.
*/
# if PRESOLVE_DEBUG

void presolve_check_reduced_costs (const CoinPostsolveMatrix *postObj)
{

  static bool warned = false ;
  static double *warned_rcosts = 0 ;
  static int allocSize = 0 ;
  static const CoinPostsolveMatrix *lastObj = 0 ;

/*
  Is the client asking for reinitialisation only?
*/
  if (postObj == 0)
  { warned = false ;
    if (warned_rcosts != 0)
    { delete[] warned_rcosts ;
      warned_rcosts = 0 ; }
    allocSize = 0 ;
    lastObj = 0 ;
    return ; }
/*
  *Should* the client have asked for reinitialisation?
*/
  int ncols0 = postObj->ncols0_ ;
  if (allocSize < ncols0 || postObj != lastObj)
  { warned = false ;
    delete[] warned_rcosts ;
    warned_rcosts = 0 ;
    allocSize = 0 ;
    lastObj = postObj ; }

  double *rcosts = postObj->rcosts_ ;

/*
  By tracking values in warned_rcosts, we can produce a single message the
  first time a value is determined to be incorrect.
*/
  if (!warned)
  { warned = true ;
    std::cout
      << "reduced cost" << std::endl ;
    warned_rcosts = new double[ncols0] ;
    CoinZeroN(warned_rcosts,ncols0) ; }

  double *colels = postObj->colels_ ;
  int *hrow = postObj->hrow_ ;
  int *mcstrt = postObj->mcstrt_ ;
  int *hincol = postObj->hincol_ ;
  CoinBigIndex *link = postObj->link_ ;

  double *clo = postObj->clo_ ;
  double *cup = postObj->cup_ ;

  double *dcost	= postObj->cost_ ;

  double *sol = postObj->sol_ ;

  char *cdone = postObj->cdone_ ;
  char *rdone = postObj->rdone_ ;

  const double ztoldj = postObj->ztoldj_ ;
  const double ztolzb = postObj->ztolzb_ ;

  double *rowduals = postObj->rowduals_ ;

  double maxmin = postObj->maxmin_ ;
  std::string strMaxmin((maxmin < 0)?"max":"min") ;
  int checkCol = -1 ;
/*
  Scan all columns, but only check the ones that are marked as having been
  postprocessed.
*/
  for (int j = 0 ; j < ncols0 ; j++)
  { if (cdone[j] == 0) continue ;
    const char *statjstr = postObj->columnStatusString(j) ;
/*
  Check the stored reduced cost for accuracy. See note above w.r.t. maxmin.
*/
    double dj = rcosts[j] ;
    double wrndj = warned_rcosts[j] ;

    { int ndx ;
      CoinBigIndex k = mcstrt[j] ;
      int len = hincol[j] ;
      double chkdj = maxmin*dcost[j] ;
      if (j == checkCol)
        std::cout
	  << "dj for " << j << " is " << dj << " - cost is " << chkdj
	  << std::endl ;
      for (ndx = 0 ; ndx < len ; ndx++)
      { int row = hrow[k] ;
	PRESOLVEASSERT(rdone[row] != 0) ;
	chkdj -= rowduals[row]*colels[k] ;
        if (j == checkCol)
	  std::cout
	    << "row " << row << " coeff " << colels[k] << " dual "
	    << rowduals[row] << " => dj " << chkdj << std::endl ;
	k = link[k] ; }
      if (fabs(dj-chkdj) > ztoldj && wrndj != dj)
      { std::cout
          << "Inacc rcost: " << j << " " << statjstr << " "
	  << strMaxmin << " have " << dj
	  << " should be " << chkdj << " err " << fabs(dj-chkdj)
	  << std::endl ; } }
/*
  Check the stored reduced cost for consistency with the variable's status.
  The cases are
    * basic: (reduced cost) == 0
    * at upper bound and not at lower bound: (reduced cost)*(maxmin) <= 0
    * at lower bound and not at upper bound: (reduced cost)*(maxmin) >= 0
    * not at either bound: any sign is correct (the variable can move either
      way) but superbasic status is sufficiently exotic that it always
      deserves a message. (There should be no superbasic variables at the
      completion of postsolve.)
  As a courtesy, show the reduced cost with the proper sign.
*/
    { double xj = sol[j] ;
      double lj = clo[j] ;
      double uj = cup[j] ;

      if (postObj->columnIsBasic(j))
      { if (fabs(dj) > ztoldj && wrndj != dj)
	{ std::cout
	    << "Bad rcost: " << j << " " << maxmin*dj
	    << " " << statjstr << " " << strMaxmin << std::endl ; } }
      else
      if (fabs(xj-uj) < ztolzb && fabs(xj-lj) > ztolzb)
      { if (dj >= ztoldj && wrndj != dj)
	{ std::cout
	    << "Bad rcost: " << j << " " << maxmin*dj
	    << " " << statjstr << " " << strMaxmin << std::endl ; } }
      else
      if (fabs(xj-lj) < ztolzb && fabs(xj-uj) > ztolzb)
      { if (dj <= -ztoldj && wrndj != dj)
	{ std::cout
	    << "Bad rcost: " << j << " " << maxmin*dj
	    << " " << statjstr << " " << strMaxmin << std::endl ; } }
      else
      if (fabs(xj-lj) > ztolzb && fabs(xj-uj) > ztolzb)
      { if (fabs(dj) > ztoldj && wrndj != dj)
        { std::cout
	    << "Superbasic rcost: " << j << " " << maxmin*dj
	    << " " << statjstr << " " << strMaxmin
	    << " lb "<< lj << " val " << xj << " ub "<< uj << std::endl ; } }
    }

    warned_rcosts[j] = rcosts[j] ; }

}

/*
  CoinPostsolveMatrix

  This routine checks the value and status of the dual variables. It
  checks that the value and status of the dual agree with the row activity.
  Errors are reported as "Bad dual"

  See presolve_check_reduced_costs for an explanation of the use of maxmin.

  Specific to CoinPostsolveMatrix due to the use of rdone. This could be fixed,
  but probably better to clone the function and specialise it for
  CoinPresolveMatrix.
*/

void presolve_check_duals (const CoinPostsolveMatrix *postObj)
{


  int nrows0 = postObj->nrows0_ ;

  double *rowduals = postObj->rowduals_ ;

  double *acts = postObj->acts_ ;
  double *rup = postObj->rup_ ;
  double *rlo = postObj->rlo_ ;

  char *rdone = postObj->rdone_ ;

  const double ztoldj = postObj->ztoldj_ ;
  const double ztolzb = postObj->ztolzb_ ;

  double maxmin = postObj->maxmin_ ;
  std::string strMaxmin((maxmin < 0)?"max":"min") ;

/*
  Scan all processed rows. The rules are as for normal reduced costs, but
  we need to remember the various flips and inversions. In summary, the correct
  situation at optimality (minimisation) is:
    * acts[i] == rup[i] ==> artificial NBLB ==> dual[i] < 0
    * acts[i] == rlo[i] ==> artificial NBUB ==> dual[i] > 0

  We can't say much about the dual for an equality. It can go either way. As a
  courtesy, show the dual with the proper sign.
*/
  for (int i = 0 ; i < nrows0 ; i++)
  { if (rdone[i] == 0) continue ;

    double ui = rup[i] ;
    double li = rlo[i] ;

    if (ui-li < 1.0e-6) continue ;

    double yi = rowduals[i] ;
    double lhsi = acts[i] ;
    const char *statistr = postObj->rowStatusString(i) ;


    if (fabs(lhsi-li) < ztolzb)
    { if (yi < -ztoldj)
      { std::cout
	  << "Bad dual: " << i << " " << maxmin*yi
	  << " " << statistr << " " << strMaxmin << std::endl ; } }
    else
    if (fabs(lhsi-ui) < ztolzb)
    { if (yi > ztoldj)
      { std::cout
	  << "Bad dual: " << i << " " << maxmin*yi
	  << " " << statistr << " " << strMaxmin << std::endl ; } }
    else
    if (li < lhsi && lhsi < ui)
    { if (fabs(yi) > ztoldj)
      { std::cout
	  << "Bad dual: " << i << " " << maxmin*yi
	  << " " << statistr << " " << strMaxmin << std::endl ; } } }
  return ; }



/*
  CoinPresolveMatrix

  This routine will check the primal (column) solution for feasibility and
  status. If there's no column solution (sol_), the routine bails out. If the
  column solution is present, all else is assumed to be present.

  chkColSol:	check colum solution (primal variables)
		0 - checks off
		1 - check for NaN/Inf
	       *2 - check for above/below column bounds
  chkRowAct:	check row solution (evaluate constraint lhs)
		0 - checks off
	       *1 - check for NaN/Inf
		2 - check for inaccuracy, above/below row bounds
  chkStatus:	check for valid status of variables
		0 - checks off
	       *1 - check status of architecturals, if colstat_ exists
	        2 - check status rows, if rowstat_ exists

  In order to check row status we need accurate row activity. Setting
  chkStatus to 2 forces chkRowAct to 2.

  CoinPrePostsolveMatrix plays games with colstat_ and rowstat_, allocating
  them as a single vector, so if colstat_ exists, rowstat_ really should
  exist. Check it anyway; this is a debug method, be robust.

  In general, the presolve transforms are not prepared to properly adjust the
  row activity (reported as `Inacc RSOL'). Postsolve transforms do better. On
  the bright side, the code seems to work just fine without maintaining row
  activity.  You probably don't want to use the level 2 checks for the row
  solution, particularly in presolve.

  With a bit of thought, the various checks could be more cleanly separated
  to require only the minimum information for each check.
*/
void presolve_check_sol (const CoinPresolveMatrix *preObj,
			 int chkColSol, int chkRowAct, int chkStatus)

{
  double *colels = preObj->colels_ ;
  int *hrow = preObj->hrow_ ;
  int *mcstrt = preObj->mcstrt_ ;
  int *hincol = preObj->hincol_ ;
  int *hinrow = preObj->hinrow_ ;

  int n	= preObj->ncols_ ;
  int m = preObj->nrows_ ;

/*
  If there's no column solution, bail out now.
*/
  if (preObj->sol_ == 0) return ;

  double *csol = preObj->sol_ ;
  double *acts = preObj->acts_ ;
  double *clo = preObj->clo_ ;
  double *cup = preObj->cup_ ;
  double *rlo = preObj->rlo_ ;
  double *rup = preObj->rup_ ;

  double tol = preObj->ztolzb_ ;

  if (chkStatus >= 2) chkRowAct = 2 ;

  double *rsol = 0 ;
  if (chkRowAct)
  { rsol = new double[m] ;
    memset(rsol,0,m*sizeof(double)) ; }

/*
  Open a loop to scan each column. For each column, do the following:
    * Update the row solution (lhs value) by adding the contribution from
      this column.
    * Check for bogus values (NaN, infinity)
    * Check for feasibility (value within column bounds)
    * Check that the status of the variable agrees with the value and with the
      lower and upper bounds. Free should have no bounds, superbasic should
      have at least one.
*/
  for (int j = 0 ; j < n ; ++j)
  { CoinBigIndex v = mcstrt[j] ;
    int colLen = hincol[j] ;
    double xj = csol[j] ;
    double lj = clo[j] ;
    double uj = cup[j] ;

    if (chkRowAct >= 1)
    { for (int u = 0 ; u < colLen ; ++u)
      { int i = hrow[v] ;
	  double aij = colels[v] ;
	  v++ ;
	  rsol[i] += aij*xj ; } }

    if (chkColSol > 0)
    { if (CoinIsnan(xj))
      { printf("NaN CSOL: %d  : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
      if (xj <= -PRESOLVE_INF || xj >= PRESOLVE_INF)
      { printf("Inf CSOL: %d  : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
      if (chkColSol > 1)
      { if (xj < lj-tol)
	{ printf("low CSOL: %d  : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
	else
	if (xj > uj+tol)
	{ printf("high CSOL: %d  : lb = %g x = %g ub = %g\n",
		 j,lj,xj,uj) ; } } }
    if (chkStatus && preObj->colstat_)
    { CoinPrePostsolveMatrix::Status statj = preObj->getColumnStatus(j) ;
      switch (statj)
      { case CoinPrePostsolveMatrix::atUpperBound:
	{ if (uj >= PRESOLVE_INF || fabs(xj-uj) > tol)
	  { printf("Bad status CSOL: %d : status atUpperBound : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::atLowerBound:
	{ if (lj <= -PRESOLVE_INF || fabs(xj-lj) > tol)
	  { printf("Bad status CSOL: %d : status atLowerBound : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::isFree:
	{ if (lj > -PRESOLVE_INF || uj < PRESOLVE_INF)
	  { printf("Bad status CSOL: %d : status isFree : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::superBasic:
	{ if (!(lj > -PRESOLVE_INF || uj < PRESOLVE_INF))
	  { printf("Bad status CSOL: %d : status superBasic : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::basic:
	{ /* Nothing to do here. */
	  break ; }
	default:
	{ printf("Bad status CSOL: %d : status unrecognized : ",j) ;
	  break ; } } } }
/*
  Now check the row solution. acts[i] is what presolve thinks we have, rsol[i]
  is what we've just calculated while scanning the columns. We need only
  check nontrivial rows (i.e., rows with length > 0). For each row,
    * Check for bogus values (NaN, infinity)
    * Check for accuracy (acts == rsol)
    * Check for feasibility (rsol within row bounds)
*/
  tol *=1.0e3;
  if (chkRowAct >= 1)
  { for (int i = 0 ; i < m ; ++i)
    { if (hinrow[i])
      { double lhsi = acts[i] ;
	double evali = rsol[i] ;
	double li = rlo[i] ;
	double ui = rup[i] ;

	if (CoinIsnan(evali) || CoinIsnan(lhsi))
	{ printf("NaN RSOL: %d  : lb = %g eval = %g (expected %g) ub = %g\n",
		 i,li,evali,lhsi,ui) ; }
	if (evali <= -PRESOLVE_INF || evali >= PRESOLVE_INF ||
	    lhsi <= -PRESOLVE_INF || lhsi >= PRESOLVE_INF)
	{ printf("Inf RSOL: %d  : lb = %g eval = %g (expected %g) ub = %g\n",
		 i,li,evali,lhsi,ui) ; }
	if (chkRowAct >= 2)
	{ if (fabs(evali-lhsi) > tol)
	  { printf("Inacc RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
		   i,li,evali,lhsi,ui) ; }
	  if (evali < li-tol || lhsi < li-tol)
	  { printf("low RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
		   i,li,evali,lhsi,ui) ; }
	  else
	  if (evali > ui+tol || lhsi > ui+tol)
	  { printf("high RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
		   i,li,evali,lhsi,ui) ; } }
	if (chkStatus >= 2 && preObj->rowstat_)
	{ CoinPrePostsolveMatrix::Status stati = preObj->getRowStatus(i) ;
	  switch (stati)
	  { case CoinPrePostsolveMatrix::atUpperBound:
	    { if (li <= -PRESOLVE_INF || fabs(lhsi-li) > tol)
	      { printf("Bad status RSOL: %d : status atUpperBound : ",i) ;
		printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ; }
	      break ; }
	    case CoinPrePostsolveMatrix::atLowerBound:
	    { if (ui >= PRESOLVE_INF || fabs(lhsi-ui) > tol)
	      { printf("Bad status RSOL: %d : status atLowerBound : ",i) ;
		printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ; }
	      break ; }
	    case CoinPrePostsolveMatrix::isFree:
	    { if (li > -PRESOLVE_INF || ui < PRESOLVE_INF)
	      { printf("Bad status RSOL: %d : status isFree : ",i) ;
		printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ; }
	      break ; }
	    case CoinPrePostsolveMatrix::superBasic:
	    { printf("Bad status RSOL: %d : status superBasic : ",i) ;
	      printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ;
	      break ; }
	    case CoinPrePostsolveMatrix::basic:
	    { /* Nothing to do here. */
	      break ; }
	    default:
	    { printf("Bad status RSOL: %d : status unrecognized : ",i) ;
	      break ; } } } } }
    delete [] rsol ; }
  return ; }

/*
  CoinPostsolveMatrix

  check_sol overload for CoinPostsolveMatrix. Parameters and functionality
  identical to check_sol immediately above, but we have to remember we're
  working with a threaded column-major representation.
*/
void presolve_check_sol (const CoinPostsolveMatrix *postObj,
			 int chkColSol, int chkRowAct, int chkStatus)

{
  double *colels = postObj->colels_ ;
  int *hrow = postObj->hrow_ ;
  int *mcstrt = postObj->mcstrt_ ;
  int *hincol = postObj->hincol_ ;
  int *link = postObj->link_ ;

  int n	= postObj->ncols_ ;
  int m = postObj->nrows_ ;

  double *csol = postObj->sol_ ;
  double *acts = postObj->acts_ ;
  double *clo = postObj->clo_ ;
  double *cup = postObj->cup_ ;
  double *rlo = postObj->rlo_ ;
  double *rup = postObj->rup_ ;

  double tol = postObj->ztolzb_ ;

  if (chkStatus >= 2) chkRowAct = 2 ;
  double *rsol = 0 ;
  if (chkRowAct >= 1)
  { rsol = new double[m] ;
    memset(rsol,0,m*sizeof(double)) ; }

/*
  Open a loop to scan each column. For each column, do the following:
    * Update the row solution (lhs value) by adding the contribution from
      this column.
    * Check for bogus values (NaN, infinity)
    * check that the status of the variable agrees with the value and with the
      lower and upper bounds. Free should have no bounds, superbasic should
      have at least one.
*/
  for (int j = 0 ; j < n ; ++j)
  { CoinBigIndex v = mcstrt[j] ;
    int colLen = hincol[j] ;
    double xj = csol[j] ;
    double lj = clo[j] ;
    double uj = cup[j] ;

    if (chkRowAct >= 1)
    { for (int u = 0 ; u < colLen ; ++u)
      { int i = hrow[v] ;
	  double aij = colels[v] ;
	  v = link[v] ;
	  rsol[i] += aij*xj ; } }
    if (chkColSol >= 1)
    { if (CoinIsnan(xj))
      { printf("NaN CSOL: %d  : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
      if (xj <= -PRESOLVE_INF || xj >= PRESOLVE_INF)
      { printf("Inf CSOL: %d  : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
      if (chkColSol >= 2)
      { if (xj < lj-tol)
	{ printf("low CSOL: %d  : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
	else
	if (xj > uj+tol)
	{ printf("high CSOL: %d  : lb = %g x = %g ub = %g\n",
		 j,lj,xj,uj) ; } } }
    if (chkStatus >= 1 && postObj->colstat_)
    { CoinPrePostsolveMatrix::Status statj = postObj->getColumnStatus(j) ;
      switch (statj)
      { case CoinPrePostsolveMatrix::atUpperBound:
	{ if (uj >= PRESOLVE_INF || fabs(xj-uj) > tol)
	  { printf("Bad status CSOL: %d : status atUpperBound : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::atLowerBound:
	{ if (lj <= -PRESOLVE_INF || fabs(xj-lj) > tol)
	  { printf("Bad status CSOL: %d : status atLowerBound : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::isFree:
	{ if (lj > -PRESOLVE_INF || uj < PRESOLVE_INF)
	  { printf("Bad status CSOL: %d : status isFree : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::superBasic:
	{ if (!(lj > -PRESOLVE_INF || uj < PRESOLVE_INF))
	  { printf("Bad status CSOL: %d : status superBasic : ",j) ;
	    printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
	  break ; }
        case CoinPrePostsolveMatrix::basic:
	{ /* Nothing to do here. */
	  break ; }
	default:
	{ printf("Bad status CSOL: %d : status unrecognized : ",j) ;
	  break ; } } } }
/*
  Now check the row solution. acts[i] is what presolve thinks we have, rsol[i]
  is what we've just calculated while scanning the columns. CoinPostsolveMatrix
  does not contain hinrow_, so we can't check for trivial rows (cf. check_sol
  for CoinPresolveMatrix).  For each row,
    * Check for bogus values (NaN, infinity)
*/
  tol *= 1.0e4;
  if (chkRowAct >= 1)
  { for (int i = 0 ; i < m ; ++i)
    { double lhsi = acts[i] ;
      double evali = rsol[i] ;
      double li = rlo[i] ;
      double ui = rup[i] ;
      
      if (CoinIsnan(evali) || CoinIsnan(lhsi))
      { printf("NaN RSOL: %d  : lb = %g eval = %g (expected %g) ub = %g\n",
	       i,li,evali,lhsi,ui) ; }
      if (evali <= -PRESOLVE_INF || evali >= PRESOLVE_INF ||
	  lhsi <= -PRESOLVE_INF || lhsi >= PRESOLVE_INF)
      { printf("Inf RSOL: %d  : lb = %g eval = %g (expected %g) ub = %g\n",
	       i,li,evali,lhsi,ui) ; }
      if (chkRowAct >= 2)
      { if (fabs(evali-lhsi) > tol)
	{ printf("Inacc RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
		 i,li,evali,lhsi,ui) ; }
        if (evali < li-tol || lhsi < li-tol)
	{ printf("low RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
		 i,li,evali,lhsi,ui) ; }
	else
	if (evali > ui+tol || lhsi > ui+tol)
	{ printf("high RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
		 i,li,evali,lhsi,ui) ; } }
      if (chkStatus >= 2 && postObj->rowstat_)
      { CoinPrePostsolveMatrix::Status stati = postObj->getRowStatus(i) ;
	switch (stati)
	{ case CoinPrePostsolveMatrix::atUpperBound:
	  { if (li <= -PRESOLVE_INF || fabs(lhsi-li) > tol)
	    { printf("Bad status RSOL: %d : status atUpperBound : ",i) ;
	      printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ; }
	    break ; }
	  case CoinPrePostsolveMatrix::atLowerBound:
	  { if (ui >= PRESOLVE_INF || fabs(lhsi-ui) > tol)
	    { printf("Bad status RSOL: %d : status atLowerBound : ",i) ;
	      printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ; }
	    break ; }
	  case CoinPrePostsolveMatrix::isFree:
	  { if (li > -PRESOLVE_INF || ui < PRESOLVE_INF)
	    { printf("Bad status RSOL: %d : status isFree : ",i) ;
	      printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ; }
	    break ; }
	  case CoinPrePostsolveMatrix::superBasic:
	  { printf("Bad status RSOL: %d : status superBasic : ",i) ;
	    printf("LB = %g lhs = %g UB = %g\n",li,lhsi,ui) ;
	    break ; }
	  case CoinPrePostsolveMatrix::basic:
	  { /* Nothing to do here. */
	    break ; }
	  default:
	  { printf("Bad status RSOL: %d : status unrecognized : ",i) ;
	    break ; } } } }
  delete [] rsol ; }
  return ; }

/*
  CoinPostsolveMatrix

  Make sure that the number of basic variables is correct.
*/
void presolve_check_nbasic (const CoinPostsolveMatrix *postObj)

{

  int ncols0 = postObj->ncols0_ ;
  int nrows0 = postObj->nrows0_ ;

  char *cdone = postObj->cdone_ ;
  char *rdone = postObj->rdone_ ;

  int nbasic = 0 ;
  int ncdone = 0;
  int nrdone = 0; 
  int ncb = 0;
  int nrb = 0;

  for (int j = 0 ; j < ncols0 ; j++)
  { 
    if (cdone[j] != 0 && postObj->columnIsBasic(j))
    { nbasic++ ;
      ncb++ ; }
    if (cdone[j])
      ncdone++ ;
  }

  for (int i = 0 ; i < nrows0 ; i++)
  {
    if (rdone[i] && postObj->rowIsBasic(i))
    { nbasic++ ;
      nrb++ ; }
    if (rdone[i])
      nrdone++ ;
  }

  if (nbasic != postObj->nrows_)
  { printf("NBASIC (ERROR): %d basic variables, should be %d; ",
	   nbasic,postObj->nrows_) ;
    printf("cdone %d, col basic %d, rdone %d, row basic %d.\n",
	   ncdone,ncb,nrdone,nrb) ;
    fflush(stdout) ; }
# if PRESOLVE_DEBUG > 1
  else
  { std::cout
      << "NBASIC: " << nbasic << " basic variables; cdone " << ncdone
      << ", col basic " << ncb << ", rdone " << nrdone << ", row basic "
      << nrb << std::endl ;
    std::cout << std::flush ;
  }
# endif
  return ; }


/*
  CoinPresolveMatrix

  Overload of presolve_check_nbasic for a CoinPresolveMatrix. There may not be
  a solution, eh?
*/
void presolve_check_nbasic (const CoinPresolveMatrix *preObj)

{

  if (preObj->sol_ == 0) return ;

  int ncols = preObj->ncols_ ;
  int nrows = preObj->nrows_ ;

  int nbasic = 0 ;
  int ncb = 0;
  int nrb = 0;

  for (int j = 0 ; j < ncols ; j++)
  { 
    if (preObj->columnIsBasic(j))
    { nbasic++ ;
      ncb++ ; }
  }

  for (int i = 0 ; i < nrows ; i++)
  {
    if (preObj->rowIsBasic(i))
    { nbasic++ ;
      nrb++ ; }
  }

  if (nbasic != nrows)
  { printf("WRONG NUMBER NBASIC:  is:  %d  should be:  %d;",
	   nbasic,nrows) ;
    printf(" cb %d, rb %d.\n",ncb,nrb);
    fflush(stdout) ; }
  return ; }

#endif
/*
  Original comment: I've forgotton what this is all about

  Looks to me like it's confirming that the columns flagged as basic indeed
  have enough coefficients between them to cover the basis. It'd be serious
  work to get this going again. Waaaaaay out of date.   -- lh, 040831 --
*/
# if 0
void check_pivots (const int *mrstrt, const int *hinrow, const int *hcol,
		   int nrows, const unsigned char *colstat,
		   const unsigned char *rowstat, int ncols)
{
  int i ;
  int nbasic = 0 ;
  int gotone = 1 ;
  int stillmore ;

  return ;

  int *bcol = new int[nrows] ;
  memset(bcol, -1, nrows*sizeof(int)) ;

  char *coldone = new char[ncols] ;
  memset(coldone, 0, ncols) ;

  while (gotone) {
    gotone = 0 ;
    stillmore = 0 ;
    for (i=0; i<nrows; i++)
      if (!postObj->rowIsBasic(i)) {
	int krs = mrstrt[i] ;
	int kre = mrstrt[i] + hinrow[i] ;
	int nb = 0 ;
	int kk ;
	for (int k=krs; k<kre; k++)
	  if (postObj->columnIsBasic(hcol[k]) && !coldone[hcol[k]]) {
	    nb++ ;
	    kk = k ;
	    if (nb > 1)
	      break ;
	  }
	if (nb == 1) {
	  PRESOLVEASSERT(bcol[i] == -1) ;
	  bcol[i] = hcol[kk] ;
	  coldone[hcol[kk]] = 1 ;
	  nbasic++ ;
	  gotone = 1 ;
	}
	else
	  stillmore = 1 ;
      }
  }
  PRESOLVEASSERT(!stillmore) ;

  for (i=0; i<nrows; i++)
    if (postObj->rowIsBasic(i)) {
      int krs = mrstrt[i] ;
      int kre = mrstrt[i] + hinrow[i] ;
      for (int k=krs; k<kre; k++)
	PRESOLVEASSERT(!postObj->columnIsBasic(hcol[k]) || coldone[hcol[k]]) ;
      nbasic++ ;
    }
  PRESOLVEASSERT(nbasic == nrows) ;
}

# endif