source: branches/devel-1/Presolve.cpp @ 35

Last change on this file since 35 was 35, checked in by forrest, 19 years ago

out CoinCopy? and CoinFill?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.8 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#define CHECK_CONSISTENCY       1
5
6#include <stdio.h>
7
8#include <assert.h>
9#include <iostream>
10
11#include "CoinHelperFunctions.hpp"
12
13#include "CoinPackedMatrix.hpp"
14#include "ClpSimplex.hpp"
15
16#include "Presolve.hpp"
17#include "PresolveMatrix.hpp"
18
19#include "PresolveEmpty.hpp"
20#include "PresolveFixed.hpp"
21#include "PresolvePsdebug.hpp"
22#include "PresolveSingleton.hpp"
23#include "PresolveDoubleton.hpp"
24#include "PresolveZeros.hpp"
25#include "PresolveSubst.hpp"
26#include "PresolveForcing.hpp"
27#include "PresolveDual.hpp"
28#include "PresolveTighten.hpp"
29#include "PresolveUseless.hpp"
30#include "PresolveDupcol.hpp"
31#include "PresolveImpliedFree.hpp"
32#include "PresolveIsolated.hpp"
33#include "ClpMessage.hpp"
34
35
36
37Presolve::Presolve() :
38  originalModel_(NULL),
39  presolvedModel_(NULL),
40  originalColumn_(NULL),
41  paction_(0),
42  ncols_(0),
43  nrows_(0),
44  nelems_(0),
45  numberPasses_(5)
46{
47}
48
49Presolve::~Presolve()
50{
51  const PresolveAction *paction = paction_;
52  while (paction) {
53    const PresolveAction *next = paction->next;
54    delete paction;
55    paction = next;
56  }
57  delete originalColumn_;
58}
59
60void Presolve::presolve(ClpSimplex& si)
61{
62  ncols_ = si.getNumCols();
63  nrows_ = si.getNumRows();
64  nelems_ = si.getNumElements();
65
66  double maxmin = si.getObjSense();
67
68  delete presolvedModel_;
69  presolvedModel_ = NULL;
70  originalModel_ = NULL;
71
72  PresolveMatrix prob(ncols_,
73                         maxmin,
74                         si,
75                         nrows_, nelems_,false);
76  prob.originalModel_ = originalModel_;
77
78  paction_ = presolve(&prob);
79
80  prob.update_model(si, nrows_, ncols_, nelems_);
81  int ncolsNow = si.getNumCols();
82  originalColumn_ = new int[ncolsNow];
83  memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
84}
85 
86
87// This is the toplevel postsolving routine.
88void Presolve::postsolve(ClpSimplex& si,
89                            unsigned char *colstat,
90                            unsigned char *rowstat)
91{
92  if (!si.isProvenOptimal()) {
93    throw CoinError("solution not optimal",
94                    "postsolve",
95                    "Presolve");
96  }
97
98  // this is the size of the original problem
99  const int ncols0  = ncols_;
100  const int nrows0  = nrows_;
101  const int nelems0 = nelems_;
102
103  double *acts = new double[nrows0];
104  double *sol  = new double[ncols0];
105
106  {
107    // put this in a block to make sure it gets deleted early
108    PostsolveMatrix prob(si,
109                            ncols0,
110                            nrows0,
111                            nelems0,
112                            si.getObjSense(),
113                            // end prepost
114
115                            sol, acts,
116                         colstat, rowstat);
117    prob.originalModel_ = originalModel_;
118
119    postsolve(prob);
120
121    // Postsolve does not use a standard matrix format, since row
122    // and col numbers change rapidly.
123    // The col rep is just a linked list.
124    int *hrow = new int[nelems0];
125    int *mcstrt = new int[ncols0+1];
126    double *colels = new double[nelems0];
127
128    {
129      const int *link = prob.link_;
130      int start = 0;
131     
132      for (int j=0; j<ncols0; j++) {
133        int k = prob.mcstrt_[j];
134        int n = prob.hincol_[j];
135
136        mcstrt[j] = start;
137        for (int i=0; i<n; i++) {
138          hrow[start] = prob.hrow_[k];
139          colels[start] = prob.colels_[k];
140          k = link[k];
141          start++;
142        }
143      }
144      mcstrt[ncols0] = start;
145    }
146   
147    si.loadProblem(ncols0, nrows0,
148                   mcstrt, hrow, colels,
149                   
150                   prob.clo_, prob.cup_,
151                   prob.cost_,
152                   prob.rlo_, prob.rup_);
153
154    delete[]hrow;
155    delete[]mcstrt;
156    delete[]colels;
157  }
158
159
160  si.setColSolution(sol);
161  delete[]acts;
162
163  delete[]sol;
164}
165/* This version of presolve returns a pointer to a new presolved
166   model.  NULL if infeasible
167*/
168ClpSimplex * 
169Presolve::presolvedModel(ClpSimplex & si,
170                         double feasibilityTolerance,
171                         bool keepIntegers,
172                         int numberPasses)
173{
174  ncols_ = si.getNumCols();
175  nrows_ = si.getNumRows();
176  nelems_ = si.getNumElements();
177  numberPasses_ = numberPasses;
178
179  double maxmin = si.getObjSense();
180  originalModel_ = &si;
181  originalColumn_ = new int[ncols_];
182
183  // Check matrix
184  if (!si.clpMatrix()->allElementsInRange(&si,1.0e-20,1.0e20))
185    return NULL;
186
187  // result is 0 - okay, 1 infeasible, -1 go round again
188  int result = -1;
189
190  while (result==-1) {
191
192    // make new copy
193    delete presolvedModel_;
194    presolvedModel_ = new ClpSimplex(si);
195   
196    // drop integer information if wanted
197    if (!keepIntegers)
198      presolvedModel_->deleteIntegerInformation();
199
200   
201    PresolveMatrix prob(ncols_,
202                        maxmin,
203                        *presolvedModel_,
204                        nrows_, nelems_,true);
205    prob.originalModel_ = originalModel_;
206
207    // move across feasibility tolerance
208    prob.feasibilityTolerance_ = feasibilityTolerance;
209
210    // Do presolve
211
212    paction_ = presolve(&prob);
213
214    result =0; // do infeasible later
215   
216    prob.update_model(*presolvedModel_, nrows_, ncols_, nelems_);
217    // copy status and solution
218    memcpy(presolvedModel_->primalColumnSolution(),
219           prob.sol_,prob.ncols_*sizeof(double));
220    memcpy(presolvedModel_->primalRowSolution(),
221      prob.acts_,prob.nrows_*sizeof(double));
222    memcpy(presolvedModel_->statusArray(),
223      prob.colstat_,prob.ncols_*sizeof(unsigned char));
224    memcpy(presolvedModel_->statusArray()+prob.ncols_,
225      prob.rowstat_,prob.nrows_*sizeof(unsigned char));
226    delete [] prob.sol_;
227    delete [] prob.acts_;
228    delete [] prob.colstat_;
229    prob.sol_=NULL;
230    prob.acts_=NULL;
231    prob.colstat_=NULL;
232
233    int ncolsNow = presolvedModel_->getNumCols();
234    memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
235    // now clean up integer variables.  This can modify original
236    int i;
237    const char * information = presolvedModel_->integerInformation();
238    if (information) {
239      int numberChanges=0;
240      double * lower0 = originalModel_->columnLower();
241      double * upper0 = originalModel_->columnUpper();
242      double * lower = presolvedModel_->columnLower();
243      double * upper = presolvedModel_->columnUpper();
244      for (i=0;i<ncolsNow;i++) {
245        if (!information[i])
246          continue;
247        int iOriginal = originalColumn_[i];
248        double lowerValue0 = lower0[iOriginal];
249        double upperValue0 = upper0[iOriginal];
250        double lowerValue = ceil(lower[i]-1.0e-5);
251        double upperValue = floor(upper[i]+1.0e-5);
252        lower[i]=lowerValue;
253        upper[i]=upperValue;
254        if (lowerValue>upperValue) {
255          numberChanges++;
256          originalModel_->messageHandler()->message(CLP_PRESOLVE_COLINFEAS,
257                                             originalModel_->messages())
258                                               <<iOriginal
259                                               <<lowerValue
260                                               <<upperValue
261                                               <<CoinMessageEol;
262          result=1;
263        } else {
264          if (lowerValue>lowerValue0+1.0e-8) {
265            lower0[iOriginal] = lowerValue;
266            numberChanges++;
267          }
268          if (upperValue<upperValue0-1.0e-8) {
269            upper0[iOriginal] = upperValue;
270            numberChanges++;
271          }
272        }         
273      }
274      if (numberChanges) {
275          originalModel_->messageHandler()->message(CLP_PRESOLVE_INTEGERMODS,
276                                             originalModel_->messages())
277                                               <<numberChanges
278                                               <<CoinMessageEol;
279        if (!result) {
280          result = -1; // round again
281        }
282      }
283    }
284  }
285  int nrowsAfter = presolvedModel_->getNumRows();
286  int ncolsAfter = presolvedModel_->getNumCols();
287  int nelsAfter = presolvedModel_->getNumElements();
288  originalModel_->messageHandler()->message(CLP_PRESOLVE_STATS,
289                                     originalModel_->messages())
290                                       <<nrowsAfter<< -(nrows_ - nrowsAfter)
291                                       << ncolsAfter<< -(ncols_ - ncolsAfter)
292                                       <<nelsAfter<< -(nelems_ - nelsAfter)
293                                       <<CoinMessageEol;
294
295  return presolvedModel_;
296}
297
298// Return pointer to presolved model
299ClpSimplex * 
300Presolve::model() const
301{
302  return presolvedModel_;
303}
304// Return pointer to original model
305ClpSimplex * 
306Presolve::originalModel() const
307{
308  return originalModel_;
309}
310void 
311Presolve::postsolve(bool updateStatus)
312{
313  if (!presolvedModel_->isProvenOptimal()) {
314    originalModel_->messageHandler()->message(CLP_PRESOLVE_NONOPTIMAL,
315                                             originalModel_->messages())
316                                               <<CoinMessageEol;
317  }
318
319  // this is the size of the original problem
320  const int ncols0  = ncols_;
321  const int nrows0  = nrows_;
322  const int nelems0 = nelems_;
323
324  // reality check
325  assert(ncols0==originalModel_->getNumCols());
326  assert(nrows0==originalModel_->getNumRows());
327
328  // this is the reduced problem
329  int ncols = presolvedModel_->getNumCols();
330  int nrows = presolvedModel_->getNumRows();
331
332  double *acts = originalModel_->primalRowSolution();
333  double *sol  = originalModel_->primalColumnSolution();
334
335  unsigned char * rowstat=NULL;
336  unsigned char * colstat = NULL;
337  if (updateStatus) {
338    unsigned char *status = originalModel_->statusArray();
339    rowstat = status + ncols0;
340    colstat = status;
341    memcpy(colstat, presolvedModel_->statusArray(), ncols);
342    memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
343  }
344
345
346  PostsolveMatrix prob(*presolvedModel_,
347                       ncols0,
348                       nrows0,
349                       nelems0,
350                       presolvedModel_->getObjSense(),
351                       // end prepost
352                       
353                       sol, acts,
354                       colstat, rowstat);
355  prob.originalModel_ = originalModel_;
356   
357  postsolve(prob);
358
359}
360
361// return pointer to original columns
362const int * 
363Presolve::originalColumns() const
364{
365  return originalColumn_;
366}
367
368// A lazy way to restrict which transformations are applied
369// during debugging.
370static int ATOI(const char *name)
371{
372  return true;
373#if     DEBUG_PRESOLVE || PRESOLVE_SUMMARY
374  if (getenv(name)) {
375    int val = atoi(getenv(name));
376    printf("%s = %d\n", name, val);
377    return (val);
378  } else {
379    if (strcmp(name,"off"))
380      return (true);
381    else
382      return (false);
383  }
384#else
385  return (true);
386#endif
387}
388
389// This is the presolve loop.
390// It is a separate virtual function so that it can be easily
391// customized by subclassing Presolve.
392const PresolveAction *Presolve::presolve(PresolveMatrix *prob)
393{
394  paction_ = 0;
395
396  prob->status_=0; // say feasible
397
398  paction_ = make_fixed(prob, paction_);
399
400#if     CHECK_CONSISTENCY
401  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
402#endif
403
404  if (!prob->status_) {
405    const bool slackd = ATOI("SLACKD");
406    //const bool forcing = ATOI("FORCING");
407    const bool doubleton = ATOI("DOUBLETON");
408    const bool forcing = ATOI("off");
409    const bool ifree = ATOI("off");
410    const bool zerocost = ATOI("off");
411    const bool dupcol = ATOI("off");
412    const bool duprow = ATOI("off");
413    const bool dual = ATOI("off");
414
415    // some things are expensive so just do once (normally)
416
417    int i;
418    // say look at all
419    for (i=0;i<nrows_;i++) 
420      prob->rowsToDo_[i]=i;
421    prob->numberRowsToDo_=nrows_;
422    for (i=0;i<ncols_;i++) 
423      prob->colsToDo_[i]=i;
424    prob->numberColsToDo_=ncols_;
425       
426
427    int iLoop;
428
429    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
430#ifdef PRESOLVE_SUMMARY
431      printf("Starting major pass %d\n",iLoop+1);
432#endif
433      const PresolveAction * const paction0 = paction_;
434      // look for substitutions with no fill
435      int fill_level=2;
436      while (1) {
437        const PresolveAction * const paction1 = paction_;
438       
439        if (slackd) {
440          bool notFinished = true;
441          while (notFinished) 
442            paction_ = slack_doubleton_action::presolve(prob, paction_,
443                                                        notFinished);
444          if (prob->status_)
445            break;
446        }
447       
448        if (doubleton) {
449          paction_ = doubleton_action::presolve(prob, paction_);
450          if (prob->status_)
451            break;
452        }
453       
454        if (zerocost) {
455          paction_ = do_tighten_action::presolve(prob, paction_);
456          if (prob->status_)
457            break;
458        }
459       
460        if (forcing) {
461          paction_ = forcing_constraint_action::presolve(prob, paction_);
462          if (prob->status_)
463            break;
464        }
465       
466        if (ifree) {
467          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
468          if (prob->status_)
469            break;
470        }
471       
472       
473#if     CHECK_CONSISTENCY
474        presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, 
475                          prob->nrows_);
476#endif
477       
478#if     DEBUG_PRESOLVE
479        presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_, 
480                          prob->ncols_);
481#endif
482#if     CHECK_CONSISTENCY
483        prob->consistent();
484#endif
485       
486#if     PRESOLVE_SUMMARY
487        prob->whichpass_++;
488#endif
489         
490        // set up for next pass
491        // later do faster if many changes i.e. memset and memcpy
492        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
493        int kcheck;
494        bool found;
495        kcheck=-1;
496        found=false;
497        for (i=0;i<prob->numberNextRowsToDo_;i++) {
498          int index = prob->nextRowsToDo_[i];
499          prob->unsetRowChanged(index);
500          prob->rowsToDo_[i] = index;
501          if (index==kcheck) {
502            printf("row %d on list after pass %d\n",kcheck,
503                   prob->whichpass_);
504            found=true;
505          }
506        }
507        if (!found&&kcheck>=0)
508          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
509        prob->numberNextRowsToDo_=0;
510        prob->numberColsToDo_ = prob->numberNextColsToDo_;
511        kcheck=-1;
512        found=false;
513        for (i=0;i<prob->numberNextColsToDo_;i++) {
514          int index = prob->nextColsToDo_[i];
515          prob->unsetColChanged(index);
516          prob->colsToDo_[i] = index;
517          if (index==kcheck) {
518            printf("col %d on list after pass %d\n",kcheck,
519                   prob->whichpass_);
520            found=true;
521          }
522        }
523        if (!found&&kcheck>=0)
524          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
525        prob->numberNextColsToDo_=0;
526        if (paction_ == paction1&&fill_level>0)
527          break;
528      }
529      // say look at all
530      int i;
531      for (i=0;i<nrows_;i++) 
532        prob->rowsToDo_[i]=i;
533      prob->numberRowsToDo_=nrows_;
534      for (i=0;i<ncols_;i++) 
535        prob->colsToDo_[i]=i;
536      prob->numberColsToDo_=ncols_;
537      // now expensive things
538      // this caused world.mps to run into numerical difficulties
539#ifdef PRESOLVE_SUMMARY
540      printf("Starting expensive\n");
541#endif
542
543      if (dual) {
544        int itry;
545        for (itry=0;itry<5;itry++) {
546          const PresolveAction * const paction2 = paction_;
547          paction_ = remove_dual_action::presolve(prob, paction_);
548          if (prob->status_)
549            break;
550          if (ifree) {
551            int fill_level=0; // switches off substitution
552            paction_ = implied_free_action::presolve(prob, paction_,fill_level);
553            if (prob->status_)
554              break;
555          }
556          if (paction_ == paction2)
557            break;
558        }
559      }
560
561      if (dupcol) {
562        paction_ = dupcol_action::presolve(prob, paction_);
563        if (prob->status_)
564          break;
565      }
566     
567      if (duprow) {
568        paction_ = duprow_action::presolve(prob, paction_);
569        if (prob->status_)
570          break;
571      }
572      if (paction_ == paction0)
573        break;
574         
575    }
576  }
577  if (!prob->status_) {
578    paction_ = drop_zero_coefficients(prob, paction_);
579
580    paction_ = drop_empty_cols_action::presolve(prob, paction_);
581    paction_ = drop_empty_rows_action::presolve(prob, paction_);
582  }
583 
584  if (prob->status_) {
585    if (prob->status_==1)
586          originalModel_->messageHandler()->message(CLP_PRESOLVE_INFEAS,
587                                             originalModel_->messages())
588                                               <<prob->feasibilityTolerance_
589                                               <<CoinMessageEol;
590    else if (prob->status_==2)
591          originalModel_->messageHandler()->message(CLP_PRESOLVE_UNBOUND,
592                                             originalModel_->messages())
593                                               <<CoinMessageEol;
594    else
595          originalModel_->messageHandler()->message(CLP_PRESOLVE_INFEASUNBOUND,
596                                             originalModel_->messages())
597                                               <<CoinMessageEol;
598    delete paction_;
599    paction_=NULL;
600  }
601  return (paction_);
602}
603
604void check_djs(PostsolveMatrix *prob);
605
606
607// We could have implemented this by having each postsolve routine
608// directly call the next one, but this may make it easier to add debugging checks.
609void Presolve::postsolve(PostsolveMatrix &prob)
610{
611  const PresolveAction *paction = paction_;
612
613  if (prob.colstat_)
614    prob.check_nbasic();
615 
616#if     DEBUG_PRESOLVE
617  check_djs(&prob);
618#endif
619 
620 
621  while (paction) {
622#if     DEBUG_PRESOLVE
623    printf("POSTSOLVING %s\n", paction->name());
624#endif
625
626    paction->postsolve(&prob);
627   
628#if     DEBUG_PRESOLVE
629    if (prob.colstat_)
630      prob.check_nbasic();
631#endif
632    paction = paction->next;
633#if     DEBUG_PRESOLVE
634    check_djs(&prob);
635#endif
636  }   
637 
638#if     0 && DEBUG_PRESOLVE
639  for (i=0; i<ncols0; i++) {
640    if (!cdone[i]) {
641      printf("!cdone[%d]\n", i);
642      abort();
643    }   
644  }
645 
646  for (i=0; i<nrows0; i++) {
647    if (!rdone[i]) {
648      printf("!rdone[%d]\n", i);
649      abort();
650    }
651  }
652 
653 
654  for (i=0; i<ncols0; i++) {
655    if (sol[i] < -1e10 || sol[i] > 1e10)
656      printf("!!!%d %g\n", i, sol[i]);
657   
658  }
659 
660 
661#endif
662 
663#if     0 && DEBUG_PRESOLVE
664  // debug check:  make sure we ended up with same original matrix
665  {
666    int identical = 1;
667   
668    for (int i=0; i<ncols0; i++) {
669      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
670      int kcs0 = &prob->mcstrt0[i];
671      int kcs = mcstrt[i];
672      int n = hincol[i];
673      for (int k=0; k<n; k++) {
674        int k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
675       
676        if (k1 == kcs+n) {
677          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
678          abort();
679        }
680       
681        if (colels[k1] != &prob->dels0[kcs0+k])
682          printf("BAD COLEL[%d %d %d]:  %g\n",
683                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
684       
685        if (kcs0+k != k1)
686          identical=0;
687      }
688    }
689    printf("identical? %d\n", identical);
690  }
691#endif
692  // put back duals
693  memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
694         nrows_*sizeof(double));
695  double maxmin = originalModel_->getObjSense();
696  if (maxmin<0.0) {
697    // swap signs
698    int i;
699    double * pi = originalModel_->dualRowSolution();
700    for (i=0;i<nrows_;i++)
701      pi[i] = -pi[i];
702  }
703  // Now check solution
704  memcpy(originalModel_->dualColumnSolution(),
705         originalModel_->objective(),ncols_*sizeof(double));
706  originalModel_->transposeTimes(-1.0,
707                                 originalModel_->dualRowSolution(),
708                                 originalModel_->dualColumnSolution());
709  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
710  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
711                        originalModel_->primalRowSolution());
712  originalModel_->checkSolution();
713  originalModel_->messageHandler()->message(CLP_PRESOLVE_POSTSOLVE,
714                                            originalModel_->messages())
715                                              <<originalModel_->objectiveValue()
716                                              <<originalModel_->sumDualInfeasibilities()
717                                              <<originalModel_->numberDualInfeasibilities()
718                                              <<originalModel_->sumPrimalInfeasibilities()
719                                              <<originalModel_->numberPrimalInfeasibilities()
720                                               <<CoinMessageEol;
721 
722}
723
724
725
726
727
728
729
730
731#if     DEBUG_PRESOLVE
732void check_djs(PostsolveMatrix *prob)
733{
734  double *colels        = prob->colels_;
735  int *hrow             = prob->hrow_;
736  int *mcstrt           = prob->mcstrt_;
737  int *hincol           = prob->hincol_;
738  int *link             = prob->link_;
739  int ncols             = prob->ncols_;
740
741  double *dcost = prob->cost_;
742
743  double *rcosts        = prob->rcosts_;
744
745  double *rowduals = prob->rowduals_;
746
747  const double maxmin   = prob->maxmin_;
748
749  char *cdone   = prob->cdone_;
750
751  double * csol = prob->sol_;
752  double * clo = prob->clo_;
753  double * cup = prob->cup_;
754  int nrows = prob->nrows_;
755  double * rlo = prob->rlo_;
756  double * rup = prob->rup_;
757  char *rdone   = prob->rdone_;
758
759  int colx;
760
761  double * rsol = new double[nrows];
762  memset(rsol,0,nrows*sizeof(double));
763
764  for (colx = 0; colx < ncols; ++colx) {
765    if (cdone[colx]) {
766      int k = mcstrt[colx];
767      int nx = hincol[colx];
768      double dj = maxmin * dcost[colx];
769      double solutionValue = csol[colx];
770      for (int i=0; i<nx; ++i) {
771        int row = hrow[k];
772        double coeff = colels[k];
773        k = link[k];
774       
775        dj -= rowduals[row] * coeff;
776        rsol[row] += solutionValue*coeff;
777      }
778      if (! (fabs(rcosts[colx] - dj) < 100*ZTOLDP))
779        printf("BAD DJ:  %d %g %g\n",
780               colx, rcosts[colx], dj);
781      if (cup[colx]-clo[colx]>1.0e-6) {
782        if (csol[colx]<clo[colx]+1.0e-6) {
783          if (dj <-1.0e-6)
784            printf("neg DJ:  %d %g  - %g %g %g\n",
785                   colx, dj, clo[colx], csol[colx], cup[colx]);
786        } else if (csol[colx]>cup[colx]-1.0e-6) {
787          if (dj > 1.0e-6)
788            printf("pos DJ:  %d %g  - %g %g %g\n",
789                   colx, dj, clo[colx], csol[colx], cup[colx]);
790        } else {
791          if (fabs(dj) >1.0e-6)
792            printf("nonzero DJ:  %d %g  - %g %g %g\n",
793                   colx, dj, clo[colx], csol[colx], cup[colx]);
794        }
795      } 
796    }
797  }
798  int rowx;
799  for (rowx = 0; rowx < nrows; ++rowx) {
800    if (rdone[rowx]) {
801      if (rup[rowx]-rlo[rowx]>1.0e-6) {
802        double dj = rowduals[rowx];
803        if (rsol[rowx]<rlo[rowx]+1.0e-6) {
804          if (dj <-1.0e-6)
805            printf("neg rDJ:  %d %g  - %g %g %g\n",
806                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
807        } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
808          if (dj > 1.0e-6)
809            printf("pos rDJ:  %d %g  - %g %g %g\n",
810                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
811        } else {
812          if (fabs(dj) >1.0e-6)
813            printf("nonzero rDJ:  %d %g  - %g %g %g\n",
814                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
815        }
816      } 
817    }
818  }
819  delete [] rsol;
820}
821#endif
822
823
Note: See TracBrowser for help on using the repository browser.