source: trunk/ClpInterior.cpp @ 275

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

Wssmp cholesky and stuff for gub

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