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

Last change on this file since 31 was 31, checked in by forrest, 17 years ago

Presolve nearly ready

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