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

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

Presolve in as option

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