source: trunk/ClpInterior.cpp @ 268

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

put back

  • 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_=2; // 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_=2; // 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_=2; // 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_=2; // 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  assert (!deltaW_);
525  deltaW_ = new double [nTotal];
526  assert (!deltaS_);
527  deltaS_ = new double [nTotal];
528  assert (!deltaT_);
529  deltaT_ = new double [nTotal];
530  assert (!zVec_);
531  zVec_ = new double [nTotal];
532  assert (!wVec_);
533  wVec_ = new double [nTotal];
534  assert (!dj_);
535  dj_ = new double [nTotal];
536  delete [] status_;
537  status_ = new unsigned char [numberRows_+numberColumns_];
538  memset(status_,0,numberRows_+numberColumns_);
539  return goodMatrix;
540}
541void
542ClpInterior::deleteWorkingData()
543{
544  int i;
545  if (optimizationDirection_!=1.0) {
546    // and modify all dual signs
547    for (i=0;i<numberColumns_;i++) 
548      reducedCost_[i] = optimizationDirection_*dj_[i];
549    for (i=0;i<numberRows_;i++) 
550      dual_[i] *= optimizationDirection_;
551  }
552  delete [] cost_;
553  cost_ = NULL;
554  delete [] solution_;
555  solution_ = NULL;
556  delete [] lower_;
557  lower_ = NULL;
558  delete [] upper_;
559  upper_ = NULL;
560  delete [] errorRegion_;
561  errorRegion_ = NULL;
562  delete [] rhsFixRegion_;
563  rhsFixRegion_ = NULL;
564  delete [] updateRegion_;
565  updateRegion_ = NULL;
566  delete [] upperSlack_;
567  upperSlack_ = NULL;
568  delete [] lowerSlack_;
569  lowerSlack_ = NULL;
570  delete [] diagonal_;
571  diagonal_ = NULL;
572  delete [] weights_;
573  weights_ = NULL;
574  delete [] deltaZ_;
575  deltaZ_ = NULL;
576  delete [] deltaW_;
577  deltaW_ = NULL;
578  delete [] deltaS_;
579  deltaS_ = NULL;
580  delete [] deltaT_;
581  deltaT_ = NULL;
582  delete [] zVec_;
583  zVec_ = NULL;
584  delete [] wVec_;
585  wVec_ = NULL;
586  delete [] dj_;
587  dj_ = NULL;
588}
589// Sanity check on input data - returns true if okay
590bool 
591ClpInterior::sanityCheck()
592{
593  // bad if empty
594  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
595    handler_->message(CLP_EMPTY_PROBLEM,messages_)
596      <<numberRows_
597      <<numberColumns_
598      <<matrix_->getNumElements()
599      <<CoinMessageEol;
600    problemStatus_=4;
601    return false;
602  }
603  int numberBad ;
604  double largestBound, smallestBound, minimumGap;
605  double smallestObj, largestObj;
606  int firstBad;
607  int modifiedBounds=0;
608  int i;
609  numberBad=0;
610  firstBad=-1;
611  minimumGap=1.0e100;
612  smallestBound=1.0e100;
613  largestBound=0.0;
614  smallestObj=1.0e100;
615  largestObj=0.0;
616  // If bounds are too close - fix
617  double fixTolerance = 1.1*primalTolerance();
618  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
619    double value;
620    value = fabs(cost_[i]);
621    if (value>1.0e50) {
622      numberBad++;
623      if (firstBad<0)
624        firstBad=i;
625    } else if (value) {
626      if (value>largestObj)
627        largestObj=value;
628      if (value<smallestObj)
629        smallestObj=value;
630    }
631    value=upper_[i]-lower_[i];
632    if (value<-primalTolerance()) {
633      numberBad++;
634      if (firstBad<0)
635        firstBad=i;
636    } else if (value<=fixTolerance) {
637      if (value) {
638        // modify
639        upper_[i] = lower_[i];
640        modifiedBounds++;
641      }
642    } else {
643      if (value<minimumGap)
644        minimumGap=value;
645    }
646    if (lower_[i]>-1.0e100&&lower_[i]) {
647      value = fabs(lower_[i]);
648      if (value>largestBound)
649        largestBound=value;
650      if (value<smallestBound)
651        smallestBound=value;
652    }
653    if (upper_[i]<1.0e100&&upper_[i]) {
654      value = fabs(upper_[i]);
655      if (value>largestBound)
656        largestBound=value;
657      if (value<smallestBound)
658        smallestBound=value;
659    }
660  }
661  if (largestBound)
662    handler_->message(CLP_RIMSTATISTICS3,messages_)
663      <<smallestBound
664      <<largestBound
665      <<minimumGap
666      <<CoinMessageEol;
667  minimumGap=1.0e100;
668  smallestBound=1.0e100;
669  largestBound=0.0;
670  for (i=0;i<numberColumns_;i++) {
671    double value;
672    value = fabs(cost_[i]);
673    if (value>1.0e50) {
674      numberBad++;
675      if (firstBad<0)
676        firstBad=i;
677    } else if (value) {
678      if (value>largestObj)
679        largestObj=value;
680      if (value<smallestObj)
681        smallestObj=value;
682    }
683    value=upper_[i]-lower_[i];
684    if (value<-primalTolerance()) {
685      numberBad++;
686      if (firstBad<0)
687        firstBad=i;
688    } else if (value<=fixTolerance) {
689      if (value) {
690        // modify
691        upper_[i] = lower_[i];
692        modifiedBounds++;
693      }
694    } else {
695      if (value<minimumGap)
696        minimumGap=value;
697    }
698    if (lower_[i]>-1.0e100&&lower_[i]) {
699      value = fabs(lower_[i]);
700      if (value>largestBound)
701        largestBound=value;
702      if (value<smallestBound)
703        smallestBound=value;
704    }
705    if (upper_[i]<1.0e100&&upper_[i]) {
706      value = fabs(upper_[i]);
707      if (value>largestBound)
708        largestBound=value;
709      if (value<smallestBound)
710        smallestBound=value;
711    }
712  }
713  char rowcol[]={'R','C'};
714  if (numberBad) {
715    handler_->message(CLP_BAD_BOUNDS,messages_)
716      <<numberBad
717      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
718      <<CoinMessageEol;
719    problemStatus_=4;
720    return false;
721  }
722  if (modifiedBounds)
723    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
724      <<modifiedBounds
725      <<CoinMessageEol;
726  handler_->message(CLP_RIMSTATISTICS1,messages_)
727    <<smallestObj
728    <<largestObj
729    <<CoinMessageEol;  if (largestBound)
730    handler_->message(CLP_RIMSTATISTICS2,messages_)
731      <<smallestBound
732      <<largestBound
733      <<minimumGap
734      <<CoinMessageEol;
735  return true;
736}
737/* Loads a problem (the constraints on the
738   rows are given by lower and upper bounds). If a pointer is 0 then the
739   following values are the default:
740   <ul>
741   <li> <code>colub</code>: all columns have upper bound infinity
742   <li> <code>collb</code>: all columns have lower bound 0
743   <li> <code>rowub</code>: all rows have upper bound infinity
744   <li> <code>rowlb</code>: all rows have lower bound -infinity
745   <li> <code>obj</code>: all variables have 0 objective coefficient
746   </ul>
747*/
748void 
749ClpInterior::loadProblem (  const ClpMatrixBase& matrix,
750                    const double* collb, const double* colub,   
751                    const double* obj,
752                    const double* rowlb, const double* rowub,
753                    const double * rowObjective)
754{
755  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
756                        rowObjective);
757}
758void 
759ClpInterior::loadProblem (  const CoinPackedMatrix& matrix,
760                    const double* collb, const double* colub,   
761                    const double* obj,
762                    const double* rowlb, const double* rowub,
763                    const double * rowObjective)
764{
765  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
766                        rowObjective);
767}
768
769/* Just like the other loadProblem() method except that the matrix is
770   given in a standard column major ordered format (without gaps). */
771void 
772ClpInterior::loadProblem (  const int numcols, const int numrows,
773                    const CoinBigIndex* start, const int* index,
774                    const double* value,
775                    const double* collb, const double* colub,   
776                    const double* obj,
777                    const double* rowlb, const double* rowub,
778                    const double * rowObjective)
779{
780  ClpModel::loadProblem(numcols, numrows, start, index, value,
781                          collb, colub, obj, rowlb, rowub,
782                          rowObjective);
783}
784void 
785ClpInterior::loadProblem (  const int numcols, const int numrows,
786                           const CoinBigIndex* start, const int* index,
787                           const double* value,const int * length,
788                           const double* collb, const double* colub,   
789                           const double* obj,
790                           const double* rowlb, const double* rowub,
791                           const double * rowObjective)
792{
793  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
794                          collb, colub, obj, rowlb, rowub,
795                          rowObjective);
796}
797// Read an mps file from the given filename
798int 
799ClpInterior::readMps(const char *filename,
800            bool keepNames,
801            bool ignoreErrors)
802{
803  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
804  return status;
805}
806#ifdef PDCO
807#include "ClpPdco.hpp"
808/* Pdco algorithm - see ClpPdco.hpp for method */
809int 
810ClpInterior::pdco()
811{
812  return ((ClpPdco *) this)->pdco();
813}
814// ** Temporary version
815int 
816ClpInterior::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
817{
818  return ((ClpPdco *) this)->pdco(stuff,options,info,outfo);
819}
820#endif
821#include "ClpPredictorCorrector.hpp"
822// Primal-Dual Predictor-Corrector barrier
823int 
824ClpInterior::primalDual()
825{ 
826  return ((ClpPredictorCorrector *) this)->solve();
827}
828
829void 
830ClpInterior::checkSolution()
831{
832  int iRow,iColumn;
833
834  objectiveValue_ = 0.0;
835  // now look at solution
836  sumPrimalInfeasibilities_=0.0;
837  sumDualInfeasibilities_=0.0;
838  double dualTolerance =  dblParam_[ClpDualTolerance];
839  double primalTolerance =  dblParam_[ClpPrimalTolerance];
840  worstComplementarity_=0.0;
841  complementarityGap_=0.0;
842
843  for (iRow=0;iRow<numberRows_;iRow++) {
844    double infeasibility=0.0;
845    double distanceUp = min(rowUpper_[iRow]-
846      rowActivity_[iRow],1.0e10);
847    double distanceDown = min(rowActivity_[iRow] -
848      rowLower_[iRow],1.0e10);
849    if (distanceUp>primalTolerance) {
850      double value = dual_[iRow];
851      // should not be negative
852      if (value<-dualTolerance) {
853        value = - value*distanceUp;
854        if (value>worstComplementarity_) 
855          worstComplementarity_=value;
856        complementarityGap_ += value;
857      }
858    }
859    if (distanceDown>primalTolerance) {
860      double value = dual_[iRow];
861      // should not be positive
862      if (value>dualTolerance) {
863        value =  value*distanceDown;
864        if (value>worstComplementarity_) 
865          worstComplementarity_=value;
866        complementarityGap_ += value;
867      }
868    }
869    if (rowActivity_[iRow]>rowUpper_[iRow]) {
870      infeasibility=rowActivity_[iRow]-rowUpper_[iRow];
871    } else if (rowActivity_[iRow]<rowLower_[iRow]) {
872      infeasibility=rowLower_[iRow]-rowActivity_[iRow];
873    }
874    if (infeasibility>primalTolerance) {
875      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
876    }
877  }
878  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
879    double infeasibility=0.0;
880    objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn];
881    double distanceUp = min(columnUpper_[iColumn]-
882      columnActivity_[iColumn],1.0e10);
883    double distanceDown = min(columnActivity_[iColumn] -
884      columnLower_[iColumn],1.0e10);
885    if (distanceUp>primalTolerance) {
886      double value = reducedCost_[iColumn];
887      // should not be negative
888      if (value<-dualTolerance) {
889        value = - value*distanceUp;
890        if (value>worstComplementarity_) 
891          worstComplementarity_=value;
892        complementarityGap_ += value;
893      }
894    }
895    if (distanceDown>primalTolerance) {
896      double value = reducedCost_[iColumn];
897      // should not be positive
898      if (value>dualTolerance) {
899        value =  value*distanceDown;
900        if (value>worstComplementarity_) 
901          worstComplementarity_=value;
902        complementarityGap_ += value;
903      }
904    }
905    if (columnActivity_[iColumn]>columnUpper_[iColumn]) {
906      infeasibility=columnActivity_[iColumn]-columnUpper_[iColumn];
907    } else if (columnActivity_[iColumn]<columnLower_[iColumn]) {
908      infeasibility=columnLower_[iColumn]-columnActivity_[iColumn];
909    }
910    if (infeasibility>primalTolerance) {
911      sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
912    }
913  }
914}
915// Set cholesky (and delete present one)
916void 
917ClpInterior::setCholesky(ClpCholeskyBase * cholesky)
918{
919  delete cholesky_;
920  cholesky_= cholesky;
921}
Note: See TracBrowser for help on using the repository browser.