source: trunk/ClpInterior.cpp @ 369

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

improving interior point code

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