// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.

#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. 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 is that the routines return void and abort when they
  find a problem.

  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
  //#define PRESOLVE_CONSISTENCY 1
  //#define PRESOLVE_DEBUG 1
#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
/*
  Check for duplicate entries in a major vector by walking the vector. For each
  coefficient, use presolve_find_row to search the remainder of the column
  for an entry with the same row 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 column 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) {
      CoinBigIndex krs = mrstrt[irow];
      CoinBigIndex kre = krs + hinrow[irow];

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

	CoinBigIndex kk = presolve_find_row1(irow, kcs, kce, hrow);
	if (kk == kce) {
	  printf("MATRIX INCONSISTENT:  can't find %s %d in %s %d\n",
		 ROW, irow, COL, jcol);
	  fflush(stdout);
	  abort();
	}
	if (testvals && colels[kk] != rowels[k]) {
	  printf("MATRIX INCONSISTENT:  value for %s %d and %s %d\n",
		 ROW, irow, COL, jcol);
	  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.
*/

void presolve_consistent(const CoinPresolveMatrix *preObj, bool testvals)
{
# if PRESOLVE_CONSISTENCY
  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") ;
# endif
}

/*
  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 PRESOLVE_CONSISTENCY
  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_) ; }
# endif

  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 PRESOLVE_CONSISTENCY
  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_) ; }
# endif

  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)
{
# if PRESOLVE_CONSISTENCY
  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();
    }
# endif
}

/*
  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 PRESOLVE_CONSISTENCY
  if (doCol)
  { links_ok(preObj->clink_,preObj->mcstrt_,
	     preObj->hincol_,preObj->ncols_) ; }
  if (doRow)
  { links_ok(preObj->rlink_,preObj->mrstrt_,
	     preObj->hinrow_,preObj->nrows_) ; }
# endif

  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)

{ 

# if PRESOLVE_CONSISTENCY
  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]) lenj-- ;

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

  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)

{ 

# if PRESOLVE_CONSISTENCY
  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) ; }
  //double *colels	= obj->colels_;
  //int *hrow		= obj->hrow_;
  CoinBigIndex *mcstrt	= obj->mcstrt_;
  int *hincol		= obj->hincol_;
  int ncols = obj->ncols_;
  for (int jcolx=0;jcolx<ncols;jcolx++) {
    CoinBigIndex k=mcstrt[jcolx];
    int i;
    for (i=0; i<hincol[jcolx]-1; ++i) {
      k = obj->link_[k];
      assert (k>=0);
    }
  }

# endif

  return ; }



/*
  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.
    * The reduced cost is calculated from scratch and compared to the
      value held in rcosts_.

  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.
*/

void presolve_check_reduced_costs (const CoinPostsolveMatrix *postObj)
{
# if PRESOLVE_DEBUG

  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 ;
/*
  Check the stored reduced cost for accuracy.
*/
    double dj = rcosts[j] ;
    double wrndj = warned_rcosts[j] ;

    { int ndx ;
      CoinBigIndex k = mcstrt[j] ;
      int len = hincol[j] ;
      double chkdj = postObj->maxmin_*dcost[j] ;
      if (j==checkCol)
        printf("dj for %d is %g - cost is %g\n",
               j,dj,chkdj);
      for (ndx = 0 ; ndx < len ; ndx++)
      { int row = hrow[k] ;
	PRESOLVEASSERT(postObj->rdone_[row] != 0) ;
	chkdj -= rowduals[row]*colels[k] ;
        if (j==checkCol)
          printf("row %d coeff %g dual %g => dj %g\n",
                 row,colels[k],rowduals[row],chkdj);

	k = link[k] ; }
      if (fabs(dj-chkdj) > ztoldj && wrndj != dj)
      { std::cout
        << j <<" "<< dj<<" " << chkdj<<" " << 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.)
*/
    { double xj = sol[j] ;
      double lj = clo[j] ;
      double uj = cup[j] ;
      const char *statjstr = postObj->columnStatusString(j) ;

      if (postObj->columnIsBasic(j))
      { if (fabs(dj) > ztoldj && wrndj != dj)
	{ std::cout
	      << j <<" "<< dj <<" "<< statjstr <<" "<< strMaxmin << std::endl ; } }
      else
      if (fabs(xj-uj) < ztolzb && fabs(xj-lj) > ztolzb)
      { if (dj >= ztoldj && wrndj != dj)
	{ std::cout
	      << j <<" "<< dj <<" "<< statjstr <<" "<< strMaxmin << std::endl ; } }
      else
      if (fabs(xj-lj) < ztolzb && fabs(xj-uj) > ztolzb)
      { if (dj <= -ztoldj && wrndj != dj)
	{ std::cout
	      << j <<" "<< dj <<" "<< statjstr <<" "<< strMaxmin << std::endl ; } }
      else
      if (fabs(xj-lj) > ztolzb && fabs(xj-uj) > ztolzb)
      { if (fabs(dj) > ztoldj && wrndj != dj)
        { std::cout
	  << j <<" "<< statjstr <<" "<< dj <<" "<< lj <<" "<< xj <<" "<< uj << std::endl ; } }
    }

    warned_rcosts[j] = rcosts[j] ; }

# endif
}

/*
  CoinPostsolveMatrix

  This routine checks the value and status of the dual variables. It
  checks that the value and status fo the dual agree with the row activity.

  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)
{
# if PRESOLVE_DEBUG


  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 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.
*/
  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
	    << i <<" "<< yi <<" "<< statistr <<" "<< strMaxmin << std::endl ; } }
    else
    if (fabs(lhsi-ui) < ztolzb)
    { if (yi > ztoldj)
      { std::cout
	    << i <<" "<< yi <<" "<< statistr <<" "<< strMaxmin << std::endl ; } }
    else
    if (li < lhsi && lhsi < ui)
    { if (fabs(yi) > ztoldj)
      { std::cout
	    << i <<" "<< yi <<" "<< statistr <<" "<< strMaxmin << std::endl ; } } }
# endif
  return ; }



/*
  This routine will check the primal (column) solution for feasibility and
  status. DON'T CALL THIS ROUTINE unless the solution is present, eh?

  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 architectural variables
		0 - checks off
	       *1 - checks on, if colstat_ exists

  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.

  To do: implement row status checks.

  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)

{
# if PRESOLVE_DEBUG
  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_ ;

  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_ ;

  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)
    { for (int u = 0 ; u < colLen ; ++u)
      { int i = hrow[v] ;
	  double aij = colels[v] ;
	  v++ ;
	  rsol[i] += aij*xj ; } }

    if (chkColSol&((1<<1)|(1<<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)
  { 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 > 1)
	{ 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) ; } } } }
    delete [] rsol ; }
# endif
  return ; }

/*
  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)

{
# if PRESOLVE_DEBUG
  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_ ;

  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 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)
    { for (int u = 0 ; u < colLen ; ++u)
      { int i = hrow[v] ;
	  double aij = colels[v] ;
	  v = link[v] ;
	  rsol[i] += aij*xj ; } }
    if (chkColSol)
    { 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 && 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)
  { 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 > 1)
      { 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) ; } } }
  delete [] rsol ; }
# endif
  return ; }

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

{
# if PRESOLVE_DEBUG

  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++ ;
    if (cdone[j])
      ncdone++;
    if (postObj->columnIsBasic(j))
      ncb++;
  }

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

  if (nbasic != postObj->nrows_)
  { printf("WRONG NUMBER NBASIC:  is:  %d  should be:  %d\n",
	   nbasic, postObj->nrows_) ;
  printf("cdone %d, cb %d, rdone %d, rb %d\n",ncdone,ncb,nrdone,nrb);
    fflush(stdout) ; }
# endif
  return ; }


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