source: trunk/ClpInterior.cpp @ 335

Last change on this file since 335 was 335, checked in by forrest, 16 years ago

solveType change

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6
7#include "CoinPragma.hpp"
8
9#include <math.h>
10
11#include "CoinHelperFunctions.hpp"
12#include "ClpInterior.hpp"
13#include "ClpMatrixBase.hpp"
14#ifdef PDCO
15#include "ClpLsqr.hpp"
16#include "ClpPdcoBase.hpp"
17#endif
18#include "CoinDenseVector.hpp"
19#include "ClpMessage.hpp"
20#include "ClpHelperFunctions.hpp"
21#include "ClpCholeskyDense.hpp"
22#include "ClpLinearObjective.hpp"
23#include <cfloat>
24
25#include <string>
26#include <stdio.h>
27#include <iostream>
28//#############################################################################
29
30ClpInterior::ClpInterior () :
31
32  ClpModel(),
33  largestPrimalError_(0.0),
34  largestDualError_(0.0),
35  sumDualInfeasibilities_(0.0),
36  sumPrimalInfeasibilities_(0.0),
37  worstComplementarity_(0.0),
38  xsize_(0.0),
39  zsize_(0.0),
40  lower_(NULL),
41  rowLowerWork_(NULL),
42  columnLowerWork_(NULL),
43  upper_(NULL),
44  rowUpperWork_(NULL),
45  columnUpperWork_(NULL),
46  cost_(NULL),
47  rhs_(NULL),
48   x_(NULL),
49   y_(NULL),
50  dj_(NULL),
51  lsqrObject_(NULL),
52  pdcoStuff_(NULL),
53  mu_(0.0),
54  objectiveNorm_(1.0e-12),
55  rhsNorm_(1.0e-12),
56  solutionNorm_(1.0e-12),
57  dualObjective_(0.0),
58  primalObjective_(0.0),
59  diagonalNorm_(1.0e-12),
60  stepLength_(0.99995),
61  linearPerturbation_(1.0e-12),
62  diagonalPerturbation_(1.0e-15),
63  targetGap_(1.0e-12),
64  projectionTolerance_(1.0e-7),
65  maximumRHSError_(0.0),
66  maximumBoundInfeasibility_(0.0),
67  maximumDualError_(0.0),
68  diagonalScaleFactor_(0.0),
69  scaleFactor_(1.0),
70  actualPrimalStep_(0.0),
71  actualDualStep_(0.0),
72  smallestInfeasibility_(0.0),
73  complementarityGap_(0.0),
74  baseObjectiveNorm_(0.0),
75  worstDirectionAccuracy_(0.0),
76  maximumRHSChange_(0.0),
77  errorRegion_(NULL),
78  rhsFixRegion_(NULL),
79  updateRegion_(NULL),
80  upperSlack_(NULL),
81  lowerSlack_(NULL),
82  diagonal_(NULL),
83  weights_(NULL),
84  solution_(NULL),
85  deltaZ_(NULL),
86  deltaW_(NULL),
87  deltaS_(NULL),
88  deltaT_(NULL),
89  zVec_(NULL),
90  wVec_(NULL),
91  cholesky_(NULL),
92  numberComplementarityPairs_(0),
93  maximumBarrierIterations_(200),
94  gonePrimalFeasible_(false),
95  goneDualFeasible_(false),
96  algorithm_(-1)
97{
98  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
99  solveType_=3; // say interior based life form
100  cholesky_ = new ClpCholeskyDense(); // put in placeholder
101}
102
103// Subproblem constructor
104ClpInterior::ClpInterior ( const ClpModel * rhs,
105                           int numberRows, const int * whichRow,
106                           int numberColumns, const int * whichColumn,
107                           bool dropNames, bool dropIntegers)
108  : ClpModel(rhs, numberRows, whichRow,
109             numberColumns,whichColumn,dropNames,dropIntegers),
110    largestPrimalError_(0.0),
111    largestDualError_(0.0),
112    sumDualInfeasibilities_(0.0),
113    sumPrimalInfeasibilities_(0.0),
114    worstComplementarity_(0.0),
115    xsize_(0.0),
116    zsize_(0.0),
117    lower_(NULL),
118    rowLowerWork_(NULL),
119    columnLowerWork_(NULL),
120    upper_(NULL),
121    rowUpperWork_(NULL),
122    columnUpperWork_(NULL),
123    cost_(NULL),
124    rhs_(NULL),
125    x_(NULL),
126    y_(NULL),
127    dj_(NULL),
128    lsqrObject_(NULL),
129    pdcoStuff_(NULL),
130    mu_(0.0),
131    objectiveNorm_(1.0e-12),
132    rhsNorm_(1.0e-12),
133    solutionNorm_(1.0e-12),
134    dualObjective_(0.0),
135    primalObjective_(0.0),
136    diagonalNorm_(1.0e-12),
137    stepLength_(0.99995),
138    linearPerturbation_(1.0e-12),
139    diagonalPerturbation_(1.0e-15),
140    targetGap_(1.0e-12),
141    projectionTolerance_(1.0e-7),
142    maximumRHSError_(0.0),
143    maximumBoundInfeasibility_(0.0),
144    maximumDualError_(0.0),
145    diagonalScaleFactor_(0.0),
146    scaleFactor_(0.0),
147    actualPrimalStep_(0.0),
148    actualDualStep_(0.0),
149    smallestInfeasibility_(0.0),
150    complementarityGap_(0.0),
151    baseObjectiveNorm_(0.0),
152    worstDirectionAccuracy_(0.0),
153    maximumRHSChange_(0.0),
154    errorRegion_(NULL),
155    rhsFixRegion_(NULL),
156    updateRegion_(NULL),
157    upperSlack_(NULL),
158    lowerSlack_(NULL),
159    diagonal_(NULL),
160    weights_(NULL),
161    solution_(NULL),
162    deltaZ_(NULL),
163    deltaW_(NULL),
164    deltaS_(NULL),
165    deltaT_(NULL),
166    zVec_(NULL),
167    wVec_(NULL),
168    cholesky_(NULL),
169    numberComplementarityPairs_(0),
170    maximumBarrierIterations_(200),
171    gonePrimalFeasible_(false),
172    goneDualFeasible_(false),
173    algorithm_(-1)
174{
175  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
176  solveType_=3; // say interior based life form
177  cholesky_= new ClpCholeskyDense();
178}
179
180//-----------------------------------------------------------------------------
181
182ClpInterior::~ClpInterior ()
183{
184  gutsOfDelete();
185}
186//#############################################################################
187/*
188   This does housekeeping
189*/
190int 
191ClpInterior::housekeeping()
192{
193  numberIterations_++;
194  return 0;
195}
196// Copy constructor.
197ClpInterior::ClpInterior(const ClpInterior &rhs) :
198  ClpModel(rhs),
199  largestPrimalError_(0.0),
200  largestDualError_(0.0),
201  sumDualInfeasibilities_(0.0),
202  sumPrimalInfeasibilities_(0.0),
203  worstComplementarity_(0.0),
204  xsize_(0.0),
205  zsize_(0.0),
206  lower_(NULL),
207  rowLowerWork_(NULL),
208  columnLowerWork_(NULL),
209  upper_(NULL),
210  rowUpperWork_(NULL),
211  columnUpperWork_(NULL),
212  cost_(NULL),
213  rhs_(NULL),
214   x_(NULL),
215   y_(NULL),
216  dj_(NULL),
217  lsqrObject_(NULL),
218  pdcoStuff_(NULL),
219  errorRegion_(NULL),
220  rhsFixRegion_(NULL),
221  updateRegion_(NULL),
222  upperSlack_(NULL),
223  lowerSlack_(NULL),
224  diagonal_(NULL),
225  weights_(NULL),
226  solution_(NULL),
227  deltaZ_(NULL),
228  deltaW_(NULL),
229  deltaS_(NULL),
230  deltaT_(NULL),
231  zVec_(NULL),
232  wVec_(NULL),
233  cholesky_(NULL)
234{
235  gutsOfDelete();
236  gutsOfCopy(rhs);
237  solveType_=3; // say interior based life form
238}
239// Copy constructor from model
240ClpInterior::ClpInterior(const ClpModel &rhs) :
241  ClpModel(rhs),
242  largestPrimalError_(0.0),
243  largestDualError_(0.0),
244  sumDualInfeasibilities_(0.0),
245  sumPrimalInfeasibilities_(0.0),
246  worstComplementarity_(0.0),
247  xsize_(0.0),
248  zsize_(0.0),
249  lower_(NULL),
250  rowLowerWork_(NULL),
251  columnLowerWork_(NULL),
252  upper_(NULL),
253  rowUpperWork_(NULL),
254  columnUpperWork_(NULL),
255  cost_(NULL),
256  rhs_(NULL),
257   x_(NULL),
258   y_(NULL),
259  dj_(NULL),
260  lsqrObject_(NULL),
261  pdcoStuff_(NULL),
262  mu_(0.0),
263  objectiveNorm_(1.0e-12),
264  rhsNorm_(1.0e-12),
265  solutionNorm_(1.0e-12),
266  dualObjective_(0.0),
267  primalObjective_(0.0),
268  diagonalNorm_(1.0e-12),
269  stepLength_(0.99995),
270  linearPerturbation_(1.0e-12),
271  diagonalPerturbation_(1.0e-15),
272  targetGap_(1.0e-12),
273  projectionTolerance_(1.0e-7),
274  maximumRHSError_(0.0),
275  maximumBoundInfeasibility_(0.0),
276  maximumDualError_(0.0),
277  diagonalScaleFactor_(0.0),
278  scaleFactor_(0.0),
279  actualPrimalStep_(0.0),
280  actualDualStep_(0.0),
281  smallestInfeasibility_(0.0),
282  complementarityGap_(0.0),
283  baseObjectiveNorm_(0.0),
284  worstDirectionAccuracy_(0.0),
285  maximumRHSChange_(0.0),
286  errorRegion_(NULL),
287  rhsFixRegion_(NULL),
288  updateRegion_(NULL),
289  upperSlack_(NULL),
290  lowerSlack_(NULL),
291  diagonal_(NULL),
292  weights_(NULL),
293  solution_(NULL),
294  deltaZ_(NULL),
295  deltaW_(NULL),
296  deltaS_(NULL),
297  deltaT_(NULL),
298  zVec_(NULL),
299  wVec_(NULL),
300  cholesky_(NULL),
301  numberComplementarityPairs_(0),
302  maximumBarrierIterations_(200),
303  gonePrimalFeasible_(false),
304  goneDualFeasible_(false),
305  algorithm_(-1)
306{
307  memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
308  solveType_=3; // say interior based life form
309  cholesky_= new ClpCholeskyDense();
310}
311// Assignment operator. This copies the data
312ClpInterior & 
313ClpInterior::operator=(const ClpInterior & rhs)
314{
315  if (this != &rhs) {
316    gutsOfDelete();
317    ClpModel::operator=(rhs);
318    gutsOfCopy(rhs);
319  }
320  return *this;
321}
322void 
323ClpInterior::gutsOfCopy(const ClpInterior & rhs)
324{
325  lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows_);
326  rowLowerWork_ = lower_+numberColumns_;
327  columnLowerWork_ = lower_;
328  upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows_);
329  rowUpperWork_ = upper_+numberColumns_;
330  columnUpperWork_ = upper_;
331  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
332  cost_ = ClpCopyOfArray(rhs.cost_,numberColumns_);
333  rhs_ = ClpCopyOfArray(rhs.rhs_,numberRows_);
334   x_ = ClpCopyOfArray(rhs.x_,numberColumns_);
335   y_ = ClpCopyOfArray(rhs.y_,numberRows_);
336  dj_ = ClpCopyOfArray(rhs.dj_,numberColumns_+numberRows_);
337#ifdef PDCO
338  lsqrObject_= new ClpLsqr(*rhs.lsqrObject_);
339  pdcoStuff_ = rhs.pdcoStuff_->clone();
340#endif
341  largestPrimalError_ = rhs.largestPrimalError_;
342  largestDualError_ = rhs.largestDualError_;
343  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
344  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
345  worstComplementarity_ = rhs.worstComplementarity_;
346  xsize_ = rhs.xsize_;
347  zsize_ = rhs.zsize_;
348  solveType_=rhs.solveType_;
349  mu_ = rhs.mu_;
350  objectiveNorm_ = rhs.objectiveNorm_;
351  rhsNorm_ = rhs.rhsNorm_;
352  solutionNorm_ = rhs.solutionNorm_;
353  dualObjective_ = rhs.dualObjective_;
354  primalObjective_ = rhs.primalObjective_;
355  diagonalNorm_ = rhs.diagonalNorm_;
356  stepLength_ = rhs.stepLength_;
357  linearPerturbation_ = rhs.linearPerturbation_;
358  diagonalPerturbation_ = rhs.diagonalPerturbation_;
359  targetGap_ = rhs.targetGap_;
360  projectionTolerance_ = rhs.projectionTolerance_;
361  maximumRHSError_ = rhs.maximumRHSError_;
362  maximumBoundInfeasibility_ = rhs.maximumBoundInfeasibility_;
363  maximumDualError_ = rhs.maximumDualError_;
364  diagonalScaleFactor_ = rhs.diagonalScaleFactor_;
365  scaleFactor_ = rhs.scaleFactor_;
366  actualPrimalStep_ = rhs.actualPrimalStep_;
367  actualDualStep_ = rhs.actualDualStep_;
368  smallestInfeasibility_ = rhs.smallestInfeasibility_;
369  complementarityGap_ = rhs.complementarityGap_;
370  baseObjectiveNorm_ = rhs.baseObjectiveNorm_;
371  worstDirectionAccuracy_ = rhs.worstDirectionAccuracy_;
372  maximumRHSChange_ = rhs.maximumRHSChange_;
373  errorRegion_ = ClpCopyOfArray(rhs.errorRegion_,numberRows_);
374  rhsFixRegion_ = ClpCopyOfArray(rhs.rhsFixRegion_,numberRows_);
375  updateRegion_ = ClpCopyOfArray(rhs.updateRegion_,numberRows_);
376  upperSlack_ = ClpCopyOfArray(rhs.upperSlack_,numberRows_+numberColumns_);
377  lowerSlack_ = ClpCopyOfArray(rhs.lowerSlack_,numberRows_+numberColumns_);
378  diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_+numberColumns_);
379  weights_ = ClpCopyOfArray(rhs.weights_,numberRows_+numberColumns_);
380  solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
381  deltaZ_ = ClpCopyOfArray(rhs.deltaZ_,numberRows_+numberColumns_);
382  deltaW_ = ClpCopyOfArray(rhs.deltaW_,numberRows_+numberColumns_);
383  deltaS_ = ClpCopyOfArray(rhs.deltaS_,numberRows_+numberColumns_);
384  deltaT_ = ClpCopyOfArray(rhs.deltaT_,numberRows_+numberColumns_);
385  zVec_ = ClpCopyOfArray(rhs.zVec_,numberRows_+numberColumns_);
386  wVec_ = ClpCopyOfArray(rhs.wVec_,numberRows_+numberColumns_);
387  cholesky_ = rhs.cholesky_->clone();
388  numberComplementarityPairs_ = rhs.numberComplementarityPairs_;
389  maximumBarrierIterations_ = rhs.maximumBarrierIterations_;
390  gonePrimalFeasible_ = rhs.gonePrimalFeasible_;
391  goneDualFeasible_ = rhs.goneDualFeasible_;
392  algorithm_ = rhs.algorithm_;
393}
394
395void 
396ClpInterior::gutsOfDelete()
397{
398  delete [] lower_;
399  lower_=NULL;
400  rowLowerWork_=NULL;
401  columnLowerWork_=NULL;
402  delete [] upper_;
403  upper_=NULL;
404  rowUpperWork_=NULL;
405  columnUpperWork_=NULL;
406  delete [] cost_;
407  cost_=NULL;
408  delete [] rhs_;
409  rhs_ = NULL;
410  delete [] x_;
411  x_ = NULL;
412  delete [] y_;
413  y_ = NULL;
414  delete [] dj_;
415  dj_ = NULL;
416#ifdef PDCO
417  delete lsqrObject_;
418  lsqrObject_ = NULL;
419  //delete pdcoStuff_;
420  pdcoStuff_=NULL;
421#endif
422  delete [] errorRegion_;
423  errorRegion_ = NULL;
424  delete [] rhsFixRegion_;
425  rhsFixRegion_ = NULL;
426  delete [] updateRegion_;
427  updateRegion_ = NULL;
428  delete [] upperSlack_;
429  upperSlack_ = NULL;
430  delete [] lowerSlack_;
431  lowerSlack_ = NULL;
432  delete [] diagonal_;
433  diagonal_ = NULL;
434  delete [] weights_;
435  weights_ = NULL;
436  delete [] solution_;
437  solution_ = NULL;
438  delete [] deltaZ_;
439  deltaZ_ = NULL;
440  delete [] deltaW_;
441  deltaW_ = NULL;
442  delete [] deltaS_;
443  deltaS_ = NULL;
444  delete [] deltaT_;
445  deltaT_ = NULL;
446  delete [] zVec_;
447  zVec_ = NULL;
448  delete [] wVec_;
449  wVec_ = NULL;
450  delete cholesky_;
451}
452bool
453ClpInterior::createWorkingData()
454{
455  bool goodMatrix=true;
456  //check matrix
457  if (!matrix_->allElementsInRange(this,1.0e-12,1.0e20)) {
458    problemStatus_=4;
459    goodMatrix= false;
460  }
461  int nTotal = numberRows_+numberColumns_;
462  delete [] solution_;
463  solution_ = new double[nTotal];
464  memcpy(solution_,columnActivity_,
465         numberColumns_*sizeof(double));
466  memcpy(solution_+numberColumns_,rowActivity_,
467         numberRows_*sizeof(double));
468  delete [] cost_;
469  cost_ = new double[nTotal];
470  int i;
471  double direction = optimizationDirection_;
472  // direction is actually scale out not scale in
473  if (direction)
474    direction = 1.0/direction;
475  const double * obj = objective();
476  for (i=0;i<numberColumns_;i++)
477    cost_[i] = direction*obj[i];
478  memset(cost_+numberColumns_,0,numberRows_*sizeof(double));
479  delete [] lower_;
480  delete [] upper_;
481  lower_ = new double[nTotal];
482  upper_ = new double[nTotal];
483  rowLowerWork_ = lower_+numberColumns_;
484  columnLowerWork_ = lower_;
485  rowUpperWork_ = upper_+numberColumns_;
486  columnUpperWork_ = upper_;
487  memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
488  memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
489  memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
490  memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
491  // clean up any mismatches on infinity
492  for (i=0;i<numberColumns_;i++) {
493    if (columnLowerWork_[i]<-1.0e30)
494      columnLowerWork_[i] = -COIN_DBL_MAX;
495    if (columnUpperWork_[i]>1.0e30)
496      columnUpperWork_[i] = COIN_DBL_MAX;
497  }
498  // clean up any mismatches on infinity
499  for (i=0;i<numberRows_;i++) {
500    if (rowLowerWork_[i]<-1.0e30)
501      rowLowerWork_[i] = -COIN_DBL_MAX;
502    if (rowUpperWork_[i]>1.0e30)
503      rowUpperWork_[i] = COIN_DBL_MAX;
504  }
505  // check rim of problem okay
506  if (!sanityCheck())
507    goodMatrix=false;
508  assert (!errorRegion_);
509  errorRegion_ = new double [numberRows_];
510  assert (!rhsFixRegion_);
511  rhsFixRegion_ = new double [numberRows_];
512  assert (!updateRegion_);
513  updateRegion_ = new double [numberRows_];
514  assert (!upperSlack_);
515  upperSlack_ = new double [nTotal];
516  assert (!lowerSlack_);
517  lowerSlack_ = new double [nTotal];
518  assert (!diagonal_);
519  diagonal_ = new double [nTotal];
520  assert (!weights_);
521  weights_ = new double [nTotal];
522  assert (!deltaZ_);
523  deltaZ_ = new double [nTotal];
524  CoinZeroN(deltaZ_,nTotal);
525  assert (!deltaW_);
526  deltaW_ = new double [nTotal];
527  CoinZeroN(deltaW_,nTotal);
528  assert (!deltaS_);
529  deltaS_ = new double [nTotal];
530  assert (!deltaT_);
531  deltaT_ = new double [nTotal];
532  assert (!zVec_);
533  zVec_ = new double [nTotal];
534  CoinZeroN(zVec_,nTotal);
535  assert (!wVec_);
536  wVec_ = new double [nTotal];
537  CoinZeroN(wVec_,nTotal);
538  assert (!dj_);
539  dj_ = new double [nTotal];
540  delete [] status_;
541  status_ = new unsigned char [numberRows_+numberColumns_];
542  memset(status_,0,numberRows_+numberColumns_);
543  return goodMatrix;
544}
545void
546ClpInterior::deleteWorkingData()
547{
548  int i;
549  if (optimizationDirection_!=1.0) {
550    // and modify all dual signs
551    for (i=0;i<numberColumns_;i++) 
552      reducedCost_[i] = optimizationDirection_*dj_[i];
553    for (i=0;i<numberRows_;i++) 
554      dual_[i] *= optimizationDirection_;
555  }
556  delete [] cost_;
557  cost_ = NULL;
558  delete [] solution_;
559  solution_ = NULL;
560  delete [] lower_;
561  lower_ = NULL;
562  delete [] upper_;
563  upper_ = NULL;
564  delete [] errorRegion_;
565  errorRegion_ = NULL;
566  delete [] rhsFixRegion_;
567  rhsFixRegion_ = NULL;
568  delete [] updateRegion_;
569  updateRegion_ = NULL;
570  delete [] upperSlack_;
571  upperSlack_ = NULL;
572  delete [] lowerSlack_;
573  lowerSlack_ = NULL;
574  delete [] diagonal_;
575  diagonal_ = NULL;
576  delete [] weights_;
577  weights_ = NULL;
578  delete [] deltaZ_;
579  deltaZ_ = NULL;
580  delete [] deltaW_;
581  deltaW_ = NULL;
582  delete [] deltaS_;
583  deltaS_ = NULL;
584  delete [] deltaT_;
585  deltaT_ = NULL;
586  delete [] zVec_;
587  zVec_ = NULL;
588  delete [] wVec_;
589  wVec_ = NULL;
590  delete [] dj_;
591  dj_ = NULL;
592}
593// Sanity check on input data - returns true if okay
594bool 
595ClpInterior::sanityCheck()
596{
597  // bad if empty
598  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
599    problemStatus_=emptyProblem();
600    return false;
601  }
602  int numberBad ;
603  double largestBound, smallestBound, minimumGap;
604  double smallestObj, largestObj;
605  int firstBad;
606  int modifiedBounds=0;
607  int i;
608  numberBad=0;
609  firstBad=-1;
610  minimumGap=1.0e100;
611  smallestBound=1.0e100;
612  largestBound=0.0;
613  smallestObj=1.0e100;
614  largestObj=0.0;
615  // If bounds are too close - fix
616  double fixTolerance = 1.1*primalTolerance();
617  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
618    double value;
619    value = fabs(cost_[i]);
620    if (value>1.0e50) {
621      numberBad++;
622      if (firstBad<0)
623        firstBad=i;
624    } else if (value) {
625      if (value>largestObj)
626        largestObj=value;
627      if (value<smallestObj)
628        smallestObj=value;
629    }
630    value=upper_[i]-lower_[i];
631    if (value<-primalTolerance()) {
632      numberBad++;
633      if (firstBad<0)
634        firstBad=i;
635    } else if (value<=fixTolerance) {
636      if (value) {
637        // modify
638        upper_[i] = lower_[i];
639        modifiedBounds++;
640      }
641    } else {
642      if (value<minimumGap)
643        minimumGap=value;
644    }
645    if (lower_[i]>-1.0e100&&lower_[i]) {
646      value = fabs(lower_[i]);
647      if (value>largestBound)
648        largestBound=value;
649      if (value<smallestBound)
650        smallestBound=value;
651    }
652    if (upper_[i]<1.0e100&&upper_[i]) {
653      value = fabs(upper_[i]);
654      if (value>largestBound)
655        largestBound=value;
656      if (value<smallestBound)
657        smallestBound=value;
658    }
659  }
660  if (largestBound)
661    handler_->message(CLP_RIMSTATISTICS3,messages_)
662      <<smallestBound
663      <<largestBound
664      <<minimumGap
665      <<CoinMessageEol;
666  minimumGap=1.0e100;
667  smallestBound=1.0e100;
668  largestBound=0.0;
669  for (i=0;i<numberColumns_;i++) {
670    double value;
671    value = fabs(cost_[i]);
672    if (value>1.0e50) {
673      numberBad++;
674      if (firstBad<0)
675        firstBad=i;
676    } else if (value) {
677      if (value>largestObj)
678        largestObj=value;
679      if (value<smallestObj)
680        smallestObj=value;
681    }
682    value=upper_[i]-lower_[i];
683    if (value<-primalTolerance()) {
684      numberBad++;
685      if (firstBad<0)
686        firstBad=i;
687    } else if (value<=fixTolerance) {
688      if (value) {
689        // modify
690        upper_[i] = lower_[i];
691        modifiedBounds++;
692      }
693    } else {
694      if (value<minimumGap)
695        minimumGap=value;
696    }
697    if (lower_[i]>-1.0e100&&lower_[i]) {
698      value = fabs(lower_[i]);
699      if (value>largestBound)
700        largestBound=value;
701      if (value<smallestBound)
702        smallestBound=value;
703    }
704    if (upper_[i]<1.0e100&&upper_[i]) {
705      value = fabs(upper_[i]);
706      if (value>largestBound)
707        largestBound=value;
708      if (value<smallestBound)
709        smallestBound=value;
710    }
711  }
712  char rowcol[]={'R','C'};
713  if (numberBad) {
714    handler_->message(CLP_BAD_BOUNDS,messages_)
715      <<numberBad
716      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
717      <<CoinMessageEol;
718    problemStatus_=4;
719    return false;
720  }
721  if (modifiedBounds)
722    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
723      <<modifiedBounds
724      <<CoinMessageEol;
725  handler_->message(CLP_RIMSTATISTICS1,messages_)
726    <<smallestObj
727    <<largestObj
728    <<CoinMessageEol;  if (largestBound)
729    handler_->message(CLP_RIMSTATISTICS2,messages_)
730      <<smallestBound
731      <<largestBound
732      <<minimumGap
733      <<CoinMessageEol;
734  return true;
735}
736/* Loads a problem (the constraints on the
737   rows are given by lower and upper bounds). If a pointer is 0 then the
738   following values are the default:
739   <ul>
740   <li> <code>colub</code>: all columns have upper bound infinity
741   <li> <code>collb</code>: all columns have lower bound 0
742   <li> <code>rowub</code>: all rows have upper bound infinity
743   <li> <code>rowlb</code>: all rows have lower bound -infinity
744   <li> <code>obj</code>: all variables have 0 objective coefficient
745   </ul>
746*/
747void 
748ClpInterior::loadProblem (  const ClpMatrixBase& matrix,
749                    const double* collb, const double* colub,   
750                    const double* obj,
751                    const double* rowlb, const double* rowub,
752                    const double * rowObjective)
753{
754  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
755                        rowObjective);
756}
757void 
758ClpInterior::loadProblem (  const CoinPackedMatrix& matrix,
759                    const double* collb, const double* colub,   
760                    const double* obj,
761                    const double* rowlb, const double* rowub,
762                    const double * rowObjective)
763{
764  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
765                        rowObjective);
766}
767
768/* Just like the other loadProblem() method except that the matrix is
769   given in a standard column major ordered format (without gaps). */
770void 
771ClpInterior::loadProblem (  const int numcols, const int numrows,
772                    const CoinBigIndex* start, const int* index,
773                    const double* value,
774                    const double* collb, const double* colub,   
775                    const double* obj,
776                    const double* rowlb, const double* rowub,
777                    const double * rowObjective)
778{
779  ClpModel::loadProblem(numcols, numrows, start, index, value,
780                          collb, colub, obj, rowlb, rowub,
781                          rowObjective);
782}
783void 
784ClpInterior::loadProblem (  const int numcols, const int numrows,
785                           const CoinBigIndex* start, const int* index,
786                           const double* value,const int * length,
787                           const double* collb, const double* colub,   
788                           const double* obj,
789                           const double* rowlb, const double* rowub,
790                           const double * rowObjective)
791{
792  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
793                          collb, colub, obj, rowlb, rowub,
794                          rowObjective);
795}
796// Read an mps file from the given filename
797int 
798ClpInterior::readMps(const char *filename,
799            bool keepNames,
800            bool ignoreErrors)
801{
802  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
803  return status;
804}
805#ifdef PDCO
806#include "ClpPdco.hpp"
807/* Pdco algorithm - see ClpPdco.hpp for method */
808int 
809ClpInterior::pdco()
810{
811  return ((ClpPdco *) this)->pdco();
812}
813// ** Temporary version
814int 
815ClpInterior::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
816{
817  return ((ClpPdco *) this)->pdco(stuff,options,info,outfo);
818}
819#endif
820#include "ClpPredictorCorrector.hpp"
821// Primal-Dual Predictor-Corrector barrier
822int 
823ClpInterior::primalDual()
824{ 
825  return ((ClpPredictorCorrector *) this)->solve();
826}
827
828void 
829ClpInterior::checkSolution()
830{
831  int iRow,iColumn;
832
833  objectiveValue_ = 0.0;
834  // now look at solution
835  sumPrimalInfeasibilities_=0.0;
836  sumDualInfeasibilities_=0.0;
837  double dualTolerance =  dblParam_[ClpDualTolerance];
838  double primalTolerance =  dblParam_[ClpPrimalTolerance];
839  worstComplementarity_=0.0;
840  complementarityGap_=0.0;
841
842  for (iRow=0;iRow<numberRows_;iRow++) {
843    double infeasibility=0.0;
844    double distanceUp = min(rowUpper_[iRow]-
845      rowActivity_[iRow],1.0e10);
846    double distanceDown = min(rowActivity_[iRow] -
847      rowLower_[iRow],1.0e10);
848    if (distanceUp>primalTolerance) {
849      double value = dual_[iRow];
850      // should not be negative
851      if (value<-dualTolerance) {
852        value = - value*distanceUp;
853        if (value>worstComplementarity_) 
854          worstComplementarity_=value;
855        complementarityGap_ += value;
856      }
857    }
858    if (distanceDown>primalTolerance) {
859      double value = dual_[iRow];
860      // should not be positive
861      if (value>dualTolerance) {
862        value =  value*distanceDown;
863        if (value>worstComplementarity_) 
864          worstComplementarity_=value;
865        complementarityGap_ += value;
866      }
867    }
868    if (rowActivity_[iRow]>rowUpper_[iRow]) {
869      infeasibility=rowActivity_[iRow]-rowUpper_[iRow];
870    } else if (rowActivity_[iRow]<rowLower_[iRow]) {
871      infeasibility=rowLower_[iRow]-rowActivity_[iRow];
872    }
873    if (infeasibility>primalTolerance) {
874      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
875    }
876  }
877  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
878    double infeasibility=0.0;
879    objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn];
880    double distanceUp = min(columnUpper_[iColumn]-
881      columnActivity_[iColumn],1.0e10);
882    double distanceDown = min(columnActivity_[iColumn] -
883      columnLower_[iColumn],1.0e10);
884    if (distanceUp>primalTolerance) {
885      double value = reducedCost_[iColumn];
886      // should not be negative
887      if (value<-dualTolerance) {
888        value = - value*distanceUp;
889        if (value>worstComplementarity_) 
890          worstComplementarity_=value;
891        complementarityGap_ += value;
892      }
893    }
894    if (distanceDown>primalTolerance) {
895      double value = reducedCost_[iColumn];
896      // should not be positive
897      if (value>dualTolerance) {
898        value =  value*distanceDown;
899        if (value>worstComplementarity_) 
900          worstComplementarity_=value;
901        complementarityGap_ += value;
902      }
903    }
904    if (columnActivity_[iColumn]>columnUpper_[iColumn]) {
905      infeasibility=columnActivity_[iColumn]-columnUpper_[iColumn];
906    } else if (columnActivity_[iColumn]<columnLower_[iColumn]) {
907      infeasibility=columnLower_[iColumn]-columnActivity_[iColumn];
908    }
909    if (infeasibility>primalTolerance) {
910      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
911    }
912  }
913}
914// Set cholesky (and delete present one)
915void 
916ClpInterior::setCholesky(ClpCholeskyBase * cholesky)
917{
918  delete cholesky_;
919  cholesky_= cholesky;
920}
Note: See TracBrowser for help on using the repository browser.