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

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

Messing about with presolve

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