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

Last change on this file since 34 was 34, checked in by forrest, 18 years ago

Allow Presolve to work with gaps

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.7 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      printf("Starting major pass %d\n",iLoop+1);
431      const PresolveAction * const paction0 = paction_;
432      // look for substitutions with no fill
433      int fill_level=2;
434      while (1) {
435        const PresolveAction * const paction1 = paction_;
436       
437        if (slackd) {
438          bool notFinished = true;
439          while (notFinished) 
440            paction_ = slack_doubleton_action::presolve(prob, paction_,
441                                                        notFinished);
442          if (prob->status_)
443            break;
444        }
445       
446        if (doubleton) {
447          paction_ = doubleton_action::presolve(prob, paction_);
448          if (prob->status_)
449            break;
450        }
451       
452        if (zerocost) {
453          paction_ = do_tighten_action::presolve(prob, paction_);
454          if (prob->status_)
455            break;
456        }
457       
458        if (forcing) {
459          paction_ = forcing_constraint_action::presolve(prob, paction_);
460          if (prob->status_)
461            break;
462        }
463       
464        if (ifree) {
465          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
466          if (prob->status_)
467            break;
468        }
469       
470       
471#if     CHECK_CONSISTENCY
472        presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, 
473                          prob->nrows_);
474#endif
475       
476#if     DEBUG_PRESOLVE
477        presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_, 
478                          prob->ncols_);
479#endif
480#if     CHECK_CONSISTENCY
481        prob->consistent();
482#endif
483       
484#if     PRESOLVE_SUMMARY
485        prob->whichpass_++;
486#endif
487         
488        // set up for next pass
489        // later do faster if many changes i.e. memset and memcpy
490        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
491        int kcheck;
492        bool found;
493        kcheck=-1;
494        found=false;
495        for (i=0;i<prob->numberNextRowsToDo_;i++) {
496          int index = prob->nextRowsToDo_[i];
497          prob->unsetRowChanged(index);
498          prob->rowsToDo_[i] = index;
499          if (index==kcheck) {
500            printf("row %d on list after pass %d\n",kcheck,
501                   prob->whichpass_);
502            found=true;
503          }
504        }
505        if (!found&&kcheck>=0)
506          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
507        prob->numberNextRowsToDo_=0;
508        prob->numberColsToDo_ = prob->numberNextColsToDo_;
509        kcheck=-1;
510        found=false;
511        for (i=0;i<prob->numberNextColsToDo_;i++) {
512          int index = prob->nextColsToDo_[i];
513          prob->unsetColChanged(index);
514          prob->colsToDo_[i] = index;
515          if (index==kcheck) {
516            printf("col %d on list after pass %d\n",kcheck,
517                   prob->whichpass_);
518            found=true;
519          }
520        }
521        if (!found&&kcheck>=0)
522          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
523        prob->numberNextColsToDo_=0;
524        if (paction_ == paction1&&fill_level>0)
525          break;
526      }
527      // say look at all
528      int i;
529      for (i=0;i<nrows_;i++) 
530        prob->rowsToDo_[i]=i;
531      prob->numberRowsToDo_=nrows_;
532      for (i=0;i<ncols_;i++) 
533        prob->colsToDo_[i]=i;
534      prob->numberColsToDo_=ncols_;
535      // now expensive things
536      // this caused world.mps to run into numerical difficulties
537      printf("Starting expensive\n");
538
539      if (dual) {
540        int itry;
541        for (itry=0;itry<5;itry++) {
542          const PresolveAction * const paction2 = paction_;
543          paction_ = remove_dual_action::presolve(prob, paction_);
544          if (prob->status_)
545            break;
546          if (ifree) {
547            int fill_level=0; // switches off substitution
548            paction_ = implied_free_action::presolve(prob, paction_,fill_level);
549            if (prob->status_)
550              break;
551          }
552          if (paction_ == paction2)
553            break;
554        }
555      }
556
557      if (dupcol) {
558        paction_ = dupcol_action::presolve(prob, paction_);
559        if (prob->status_)
560          break;
561      }
562     
563      if (duprow) {
564        paction_ = duprow_action::presolve(prob, paction_);
565        if (prob->status_)
566          break;
567      }
568      if (paction_ == paction0)
569        break;
570         
571    }
572  }
573  if (!prob->status_) {
574    paction_ = drop_zero_coefficients(prob, paction_);
575
576    paction_ = drop_empty_cols_action::presolve(prob, paction_);
577    paction_ = drop_empty_rows_action::presolve(prob, paction_);
578  }
579 
580  if (prob->status_) {
581    if (prob->status_==1)
582          originalModel_->messageHandler()->message(CLP_PRESOLVE_INFEAS,
583                                             originalModel_->messages())
584                                               <<prob->feasibilityTolerance_
585                                               <<CoinMessageEol;
586    else if (prob->status_==2)
587          originalModel_->messageHandler()->message(CLP_PRESOLVE_UNBOUND,
588                                             originalModel_->messages())
589                                               <<CoinMessageEol;
590    else
591          originalModel_->messageHandler()->message(CLP_PRESOLVE_INFEASUNBOUND,
592                                             originalModel_->messages())
593                                               <<CoinMessageEol;
594    delete paction_;
595    paction_=NULL;
596  }
597  return (paction_);
598}
599
600void check_djs(PostsolveMatrix *prob);
601
602
603// We could have implemented this by having each postsolve routine
604// directly call the next one, but this may make it easier to add debugging checks.
605void Presolve::postsolve(PostsolveMatrix &prob)
606{
607  const PresolveAction *paction = paction_;
608
609  if (prob.colstat_)
610    prob.check_nbasic();
611 
612#if     DEBUG_PRESOLVE
613  check_djs(&prob);
614#endif
615 
616 
617  while (paction) {
618#if     DEBUG_PRESOLVE
619    printf("POSTSOLVING %s\n", paction->name());
620#endif
621
622    paction->postsolve(&prob);
623   
624#if     DEBUG_PRESOLVE
625    if (prob.colstat_)
626      prob.check_nbasic();
627#endif
628    paction = paction->next;
629#if     DEBUG_PRESOLVE
630    check_djs(&prob);
631#endif
632  }   
633 
634#if     0 && DEBUG_PRESOLVE
635  for (i=0; i<ncols0; i++) {
636    if (!cdone[i]) {
637      printf("!cdone[%d]\n", i);
638      abort();
639    }   
640  }
641 
642  for (i=0; i<nrows0; i++) {
643    if (!rdone[i]) {
644      printf("!rdone[%d]\n", i);
645      abort();
646    }
647  }
648 
649 
650  for (i=0; i<ncols0; i++) {
651    if (sol[i] < -1e10 || sol[i] > 1e10)
652      printf("!!!%d %g\n", i, sol[i]);
653   
654  }
655 
656 
657#endif
658 
659#if     0 && DEBUG_PRESOLVE
660  // debug check:  make sure we ended up with same original matrix
661  {
662    int identical = 1;
663   
664    for (int i=0; i<ncols0; i++) {
665      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
666      int kcs0 = &prob->mcstrt0[i];
667      int kcs = mcstrt[i];
668      int n = hincol[i];
669      for (int k=0; k<n; k++) {
670        int k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
671       
672        if (k1 == kcs+n) {
673          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
674          abort();
675        }
676       
677        if (colels[k1] != &prob->dels0[kcs0+k])
678          printf("BAD COLEL[%d %d %d]:  %g\n",
679                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
680       
681        if (kcs0+k != k1)
682          identical=0;
683      }
684    }
685    printf("identical? %d\n", identical);
686  }
687#endif
688  // put back duals
689  memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
690         nrows_*sizeof(double));
691  double maxmin = originalModel_->getObjSense();
692  if (maxmin<0.0) {
693    // swap signs
694    int i;
695    double * pi = originalModel_->dualRowSolution();
696    for (i=0;i<nrows_;i++)
697      pi[i] = -pi[i];
698  }
699  // Now check solution
700  memcpy(originalModel_->dualColumnSolution(),
701         originalModel_->objective(),ncols_*sizeof(double));
702  originalModel_->transposeTimes(-1.0,
703                                 originalModel_->dualRowSolution(),
704                                 originalModel_->dualColumnSolution());
705  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
706  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
707                        originalModel_->primalRowSolution());
708  originalModel_->checkSolution();
709  originalModel_->messageHandler()->message(CLP_PRESOLVE_POSTSOLVE,
710                                            originalModel_->messages())
711                                              <<originalModel_->objectiveValue()
712                                              <<originalModel_->sumDualInfeasibilities()
713                                              <<originalModel_->numberDualInfeasibilities()
714                                              <<originalModel_->sumPrimalInfeasibilities()
715                                              <<originalModel_->numberPrimalInfeasibilities()
716                                               <<CoinMessageEol;
717 
718}
719
720
721
722
723
724
725
726
727#if     DEBUG_PRESOLVE
728void check_djs(PostsolveMatrix *prob)
729{
730  double *colels        = prob->colels_;
731  int *hrow             = prob->hrow_;
732  int *mcstrt           = prob->mcstrt_;
733  int *hincol           = prob->hincol_;
734  int *link             = prob->link_;
735  int ncols             = prob->ncols_;
736
737  double *dcost = prob->cost_;
738
739  double *rcosts        = prob->rcosts_;
740
741  double *rowduals = prob->rowduals_;
742
743  const double maxmin   = prob->maxmin_;
744
745  char *cdone   = prob->cdone_;
746
747  double * csol = prob->sol_;
748  double * clo = prob->clo_;
749  double * cup = prob->cup_;
750  int nrows = prob->nrows_;
751  double * rlo = prob->rlo_;
752  double * rup = prob->rup_;
753  char *rdone   = prob->rdone_;
754
755  int colx;
756
757  double * rsol = new double[nrows];
758  memset(rsol,0,nrows*sizeof(double));
759
760  for (colx = 0; colx < ncols; ++colx) {
761    if (cdone[colx]) {
762      int k = mcstrt[colx];
763      int nx = hincol[colx];
764      double dj = maxmin * dcost[colx];
765      double solutionValue = csol[colx];
766      for (int i=0; i<nx; ++i) {
767        int row = hrow[k];
768        double coeff = colels[k];
769        k = link[k];
770       
771        dj -= rowduals[row] * coeff;
772        rsol[row] += solutionValue*coeff;
773      }
774      if (! (fabs(rcosts[colx] - dj) < 100*ZTOLDP))
775        printf("BAD DJ:  %d %g %g\n",
776               colx, rcosts[colx], dj);
777      if (cup[colx]-clo[colx]>1.0e-6) {
778        if (csol[colx]<clo[colx]+1.0e-6) {
779          if (dj <-1.0e-6)
780            printf("neg DJ:  %d %g  - %g %g %g\n",
781                   colx, dj, clo[colx], csol[colx], cup[colx]);
782        } else if (csol[colx]>cup[colx]-1.0e-6) {
783          if (dj > 1.0e-6)
784            printf("pos DJ:  %d %g  - %g %g %g\n",
785                   colx, dj, clo[colx], csol[colx], cup[colx]);
786        } else {
787          if (fabs(dj) >1.0e-6)
788            printf("nonzero DJ:  %d %g  - %g %g %g\n",
789                   colx, dj, clo[colx], csol[colx], cup[colx]);
790        }
791      } 
792    }
793  }
794  int rowx;
795  for (rowx = 0; rowx < nrows; ++rowx) {
796    if (rdone[rowx]) {
797      if (rup[rowx]-rlo[rowx]>1.0e-6) {
798        double dj = rowduals[rowx];
799        if (rsol[rowx]<rlo[rowx]+1.0e-6) {
800          if (dj <-1.0e-6)
801            printf("neg rDJ:  %d %g  - %g %g %g\n",
802                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
803        } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
804          if (dj > 1.0e-6)
805            printf("pos rDJ:  %d %g  - %g %g %g\n",
806                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
807        } else {
808          if (fabs(dj) >1.0e-6)
809            printf("nonzero rDJ:  %d %g  - %g %g %g\n",
810                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
811        }
812      } 
813    }
814  }
815  delete [] rsol;
816}
817#endif
818
819
Note: See TracBrowser for help on using the repository browser.