source: trunk/Clp/src/ClpInterior.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 37.3 KB
Line 
1/* $Id: ClpInterior.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7
8#include <math.h>
9
10#include "CoinHelperFunctions.hpp"
11#include "ClpInterior.hpp"
12#include "ClpMatrixBase.hpp"
13#include "ClpLsqr.hpp"
14#include "ClpPdcoBase.hpp"
15#include "CoinDenseVector.hpp"
16#include "ClpMessage.hpp"
17#include "ClpHelperFunctions.hpp"
18#include "ClpCholeskyDense.hpp"
19#include "ClpLinearObjective.hpp"
20#include "ClpQuadraticObjective.hpp"
21#include <cfloat>
22
23#include <string>
24#include <stdio.h>
25#include <iostream>
26//#############################################################################
27
28ClpInterior::ClpInterior()
29  :
30
31  ClpModel()
32  , largestPrimalError_(0.0)
33  , largestDualError_(0.0)
34  , sumDualInfeasibilities_(0.0)
35  , sumPrimalInfeasibilities_(0.0)
36  , worstComplementarity_(0.0)
37  , xsize_(0.0)
38  , zsize_(0.0)
39  , lower_(NULL)
40  , rowLowerWork_(NULL)
41  , columnLowerWork_(NULL)
42  , upper_(NULL)
43  , rowUpperWork_(NULL)
44  , columnUpperWork_(NULL)
45  , cost_(NULL)
46  , rhs_(NULL)
47  , x_(NULL)
48  , y_(NULL)
49  , dj_(NULL)
50  , lsqrObject_(NULL)
51  , pdcoStuff_(NULL)
52  , mu_(0.0)
53  , objectiveNorm_(1.0e-12)
54  , rhsNorm_(1.0e-12)
55  , solutionNorm_(1.0e-12)
56  , dualObjective_(0.0)
57  , primalObjective_(0.0)
58  , diagonalNorm_(1.0e-12)
59  , stepLength_(0.995)
60  , linearPerturbation_(1.0e-12)
61  , diagonalPerturbation_(1.0e-15)
62  , gamma_(0.0)
63  , delta_(0)
64  , targetGap_(1.0e-12)
65  , projectionTolerance_(1.0e-7)
66  , maximumRHSError_(0.0)
67  , maximumBoundInfeasibility_(0.0)
68  , maximumDualError_(0.0)
69  , diagonalScaleFactor_(0.0)
70  , scaleFactor_(1.0)
71  , actualPrimalStep_(0.0)
72  , actualDualStep_(0.0)
73  , smallestInfeasibility_(0.0)
74  , complementarityGap_(0.0)
75  , baseObjectiveNorm_(0.0)
76  , worstDirectionAccuracy_(0.0)
77  , maximumRHSChange_(0.0)
78  , errorRegion_(NULL)
79  , rhsFixRegion_(NULL)
80  , upperSlack_(NULL)
81  , lowerSlack_(NULL)
82  , diagonal_(NULL)
83  , solution_(NULL)
84  , workArray_(NULL)
85  , deltaX_(NULL)
86  , deltaY_(NULL)
87  , deltaZ_(NULL)
88  , deltaW_(NULL)
89  , deltaSU_(NULL)
90  , deltaSL_(NULL)
91  , primalR_(NULL)
92  , dualR_(NULL)
93  , rhsB_(NULL)
94  , rhsU_(NULL)
95  , rhsL_(NULL)
96  , rhsZ_(NULL)
97  , rhsW_(NULL)
98  , rhsC_(NULL)
99  , zVec_(NULL)
100  , wVec_(NULL)
101  , cholesky_(NULL)
102  , numberComplementarityPairs_(0)
103  , numberComplementarityItems_(0)
104  , maximumBarrierIterations_(200)
105  , gonePrimalFeasible_(false)
106  , goneDualFeasible_(false)
107  , algorithm_(-1)
108{
109  memset(historyInfeasibility_, 0, LENGTH_HISTORY * sizeof(CoinWorkDouble));
110  solveType_ = 3; // say interior based life form
111  cholesky_ = new ClpCholeskyDense(); // put in placeholder
112}
113
114// Subproblem constructor
115ClpInterior::ClpInterior(const ClpModel *rhs,
116  int numberRows, const int *whichRow,
117  int numberColumns, const int *whichColumn,
118  bool dropNames, bool dropIntegers)
119  : ClpModel(rhs, numberRows, whichRow,
120      numberColumns, whichColumn, dropNames, dropIntegers)
121  , largestPrimalError_(0.0)
122  , largestDualError_(0.0)
123  , sumDualInfeasibilities_(0.0)
124  , sumPrimalInfeasibilities_(0.0)
125  , worstComplementarity_(0.0)
126  , xsize_(0.0)
127  , zsize_(0.0)
128  , lower_(NULL)
129  , rowLowerWork_(NULL)
130  , columnLowerWork_(NULL)
131  , upper_(NULL)
132  , rowUpperWork_(NULL)
133  , columnUpperWork_(NULL)
134  , cost_(NULL)
135  , rhs_(NULL)
136  , x_(NULL)
137  , y_(NULL)
138  , dj_(NULL)
139  , lsqrObject_(NULL)
140  , pdcoStuff_(NULL)
141  , mu_(0.0)
142  , objectiveNorm_(1.0e-12)
143  , rhsNorm_(1.0e-12)
144  , solutionNorm_(1.0e-12)
145  , dualObjective_(0.0)
146  , primalObjective_(0.0)
147  , diagonalNorm_(1.0e-12)
148  , stepLength_(0.99995)
149  , linearPerturbation_(1.0e-12)
150  , diagonalPerturbation_(1.0e-15)
151  , gamma_(0.0)
152  , delta_(0)
153  , targetGap_(1.0e-12)
154  , projectionTolerance_(1.0e-7)
155  , maximumRHSError_(0.0)
156  , maximumBoundInfeasibility_(0.0)
157  , maximumDualError_(0.0)
158  , diagonalScaleFactor_(0.0)
159  , scaleFactor_(0.0)
160  , actualPrimalStep_(0.0)
161  , actualDualStep_(0.0)
162  , smallestInfeasibility_(0.0)
163  , complementarityGap_(0.0)
164  , baseObjectiveNorm_(0.0)
165  , worstDirectionAccuracy_(0.0)
166  , maximumRHSChange_(0.0)
167  , errorRegion_(NULL)
168  , rhsFixRegion_(NULL)
169  , upperSlack_(NULL)
170  , lowerSlack_(NULL)
171  , diagonal_(NULL)
172  , solution_(NULL)
173  , workArray_(NULL)
174  , deltaX_(NULL)
175  , deltaY_(NULL)
176  , deltaZ_(NULL)
177  , deltaW_(NULL)
178  , deltaSU_(NULL)
179  , deltaSL_(NULL)
180  , primalR_(NULL)
181  , dualR_(NULL)
182  , rhsB_(NULL)
183  , rhsU_(NULL)
184  , rhsL_(NULL)
185  , rhsZ_(NULL)
186  , rhsW_(NULL)
187  , rhsC_(NULL)
188  , zVec_(NULL)
189  , wVec_(NULL)
190  , cholesky_(NULL)
191  , numberComplementarityPairs_(0)
192  , numberComplementarityItems_(0)
193  , maximumBarrierIterations_(200)
194  , gonePrimalFeasible_(false)
195  , goneDualFeasible_(false)
196  , algorithm_(-1)
197{
198  memset(historyInfeasibility_, 0, LENGTH_HISTORY * sizeof(CoinWorkDouble));
199  solveType_ = 3; // say interior based life form
200  cholesky_ = new ClpCholeskyDense();
201}
202
203//-----------------------------------------------------------------------------
204
205ClpInterior::~ClpInterior()
206{
207  gutsOfDelete();
208}
209//#############################################################################
210/*
211   This does housekeeping
212*/
213int ClpInterior::housekeeping()
214{
215  numberIterations_++;
216  return 0;
217}
218// Copy constructor.
219ClpInterior::ClpInterior(const ClpInterior &rhs)
220  : ClpModel(rhs)
221  , largestPrimalError_(0.0)
222  , largestDualError_(0.0)
223  , sumDualInfeasibilities_(0.0)
224  , sumPrimalInfeasibilities_(0.0)
225  , worstComplementarity_(0.0)
226  , xsize_(0.0)
227  , zsize_(0.0)
228  , lower_(NULL)
229  , rowLowerWork_(NULL)
230  , columnLowerWork_(NULL)
231  , upper_(NULL)
232  , rowUpperWork_(NULL)
233  , columnUpperWork_(NULL)
234  , cost_(NULL)
235  , rhs_(NULL)
236  , x_(NULL)
237  , y_(NULL)
238  , dj_(NULL)
239  , lsqrObject_(NULL)
240  , pdcoStuff_(NULL)
241  , errorRegion_(NULL)
242  , rhsFixRegion_(NULL)
243  , upperSlack_(NULL)
244  , lowerSlack_(NULL)
245  , diagonal_(NULL)
246  , solution_(NULL)
247  , workArray_(NULL)
248  , deltaX_(NULL)
249  , deltaY_(NULL)
250  , deltaZ_(NULL)
251  , deltaW_(NULL)
252  , deltaSU_(NULL)
253  , deltaSL_(NULL)
254  , primalR_(NULL)
255  , dualR_(NULL)
256  , rhsB_(NULL)
257  , rhsU_(NULL)
258  , rhsL_(NULL)
259  , rhsZ_(NULL)
260  , rhsW_(NULL)
261  , rhsC_(NULL)
262  , zVec_(NULL)
263  , wVec_(NULL)
264  , cholesky_(NULL)
265{
266  gutsOfDelete();
267  gutsOfCopy(rhs);
268  solveType_ = 3; // say interior based life form
269}
270// Copy constructor from model
271ClpInterior::ClpInterior(const ClpModel &rhs)
272  : ClpModel(rhs)
273  , largestPrimalError_(0.0)
274  , largestDualError_(0.0)
275  , sumDualInfeasibilities_(0.0)
276  , sumPrimalInfeasibilities_(0.0)
277  , worstComplementarity_(0.0)
278  , xsize_(0.0)
279  , zsize_(0.0)
280  , lower_(NULL)
281  , rowLowerWork_(NULL)
282  , columnLowerWork_(NULL)
283  , upper_(NULL)
284  , rowUpperWork_(NULL)
285  , columnUpperWork_(NULL)
286  , cost_(NULL)
287  , rhs_(NULL)
288  , x_(NULL)
289  , y_(NULL)
290  , dj_(NULL)
291  , lsqrObject_(NULL)
292  , pdcoStuff_(NULL)
293  , mu_(0.0)
294  , objectiveNorm_(1.0e-12)
295  , rhsNorm_(1.0e-12)
296  , solutionNorm_(1.0e-12)
297  , dualObjective_(0.0)
298  , primalObjective_(0.0)
299  , diagonalNorm_(1.0e-12)
300  , stepLength_(0.99995)
301  , linearPerturbation_(1.0e-12)
302  , diagonalPerturbation_(1.0e-15)
303  , gamma_(0.0)
304  , delta_(0)
305  , targetGap_(1.0e-12)
306  , projectionTolerance_(1.0e-7)
307  , maximumRHSError_(0.0)
308  , maximumBoundInfeasibility_(0.0)
309  , maximumDualError_(0.0)
310  , diagonalScaleFactor_(0.0)
311  , scaleFactor_(0.0)
312  , actualPrimalStep_(0.0)
313  , actualDualStep_(0.0)
314  , smallestInfeasibility_(0.0)
315  , complementarityGap_(0.0)
316  , baseObjectiveNorm_(0.0)
317  , worstDirectionAccuracy_(0.0)
318  , maximumRHSChange_(0.0)
319  , errorRegion_(NULL)
320  , rhsFixRegion_(NULL)
321  , upperSlack_(NULL)
322  , lowerSlack_(NULL)
323  , diagonal_(NULL)
324  , solution_(NULL)
325  , workArray_(NULL)
326  , deltaX_(NULL)
327  , deltaY_(NULL)
328  , deltaZ_(NULL)
329  , deltaW_(NULL)
330  , deltaSU_(NULL)
331  , deltaSL_(NULL)
332  , primalR_(NULL)
333  , dualR_(NULL)
334  , rhsB_(NULL)
335  , rhsU_(NULL)
336  , rhsL_(NULL)
337  , rhsZ_(NULL)
338  , rhsW_(NULL)
339  , rhsC_(NULL)
340  , zVec_(NULL)
341  , wVec_(NULL)
342  , cholesky_(NULL)
343  , numberComplementarityPairs_(0)
344  , numberComplementarityItems_(0)
345  , maximumBarrierIterations_(200)
346  , gonePrimalFeasible_(false)
347  , goneDualFeasible_(false)
348  , algorithm_(-1)
349{
350  memset(historyInfeasibility_, 0, LENGTH_HISTORY * sizeof(CoinWorkDouble));
351  solveType_ = 3; // say interior based life form
352  cholesky_ = new ClpCholeskyDense();
353}
354// Assignment operator. This copies the data
355ClpInterior &
356ClpInterior::operator=(const ClpInterior &rhs)
357{
358  if (this != &rhs) {
359    gutsOfDelete();
360    ClpModel::operator=(rhs);
361    gutsOfCopy(rhs);
362  }
363  return *this;
364}
365void ClpInterior::gutsOfCopy(const ClpInterior &rhs)
366{
367  lower_ = ClpCopyOfArray(rhs.lower_, numberColumns_ + numberRows_);
368  rowLowerWork_ = lower_ + numberColumns_;
369  columnLowerWork_ = lower_;
370  upper_ = ClpCopyOfArray(rhs.upper_, numberColumns_ + numberRows_);
371  rowUpperWork_ = upper_ + numberColumns_;
372  columnUpperWork_ = upper_;
373  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
374  cost_ = ClpCopyOfArray(rhs.cost_, numberColumns_);
375  rhs_ = ClpCopyOfArray(rhs.rhs_, numberRows_);
376  x_ = ClpCopyOfArray(rhs.x_, numberColumns_);
377  y_ = ClpCopyOfArray(rhs.y_, numberRows_);
378  dj_ = ClpCopyOfArray(rhs.dj_, numberColumns_ + numberRows_);
379  lsqrObject_ = rhs.lsqrObject_ != NULL ? new ClpLsqr(*rhs.lsqrObject_) : NULL;
380  pdcoStuff_ = rhs.pdcoStuff_ != NULL ? rhs.pdcoStuff_->clone() : NULL;
381  largestPrimalError_ = rhs.largestPrimalError_;
382  largestDualError_ = rhs.largestDualError_;
383  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
384  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
385  worstComplementarity_ = rhs.worstComplementarity_;
386  xsize_ = rhs.xsize_;
387  zsize_ = rhs.zsize_;
388  solveType_ = rhs.solveType_;
389  mu_ = rhs.mu_;
390  objectiveNorm_ = rhs.objectiveNorm_;
391  rhsNorm_ = rhs.rhsNorm_;
392  solutionNorm_ = rhs.solutionNorm_;
393  dualObjective_ = rhs.dualObjective_;
394  primalObjective_ = rhs.primalObjective_;
395  diagonalNorm_ = rhs.diagonalNorm_;
396  stepLength_ = rhs.stepLength_;
397  linearPerturbation_ = rhs.linearPerturbation_;
398  diagonalPerturbation_ = rhs.diagonalPerturbation_;
399  gamma_ = rhs.gamma_;
400  delta_ = rhs.delta_;
401  targetGap_ = rhs.targetGap_;
402  projectionTolerance_ = rhs.projectionTolerance_;
403  maximumRHSError_ = rhs.maximumRHSError_;
404  maximumBoundInfeasibility_ = rhs.maximumBoundInfeasibility_;
405  maximumDualError_ = rhs.maximumDualError_;
406  diagonalScaleFactor_ = rhs.diagonalScaleFactor_;
407  scaleFactor_ = rhs.scaleFactor_;
408  actualPrimalStep_ = rhs.actualPrimalStep_;
409  actualDualStep_ = rhs.actualDualStep_;
410  smallestInfeasibility_ = rhs.smallestInfeasibility_;
411  complementarityGap_ = rhs.complementarityGap_;
412  baseObjectiveNorm_ = rhs.baseObjectiveNorm_;
413  worstDirectionAccuracy_ = rhs.worstDirectionAccuracy_;
414  maximumRHSChange_ = rhs.maximumRHSChange_;
415  errorRegion_ = ClpCopyOfArray(rhs.errorRegion_, numberRows_);
416  rhsFixRegion_ = ClpCopyOfArray(rhs.rhsFixRegion_, numberRows_);
417  deltaY_ = ClpCopyOfArray(rhs.deltaY_, numberRows_);
418  upperSlack_ = ClpCopyOfArray(rhs.upperSlack_, numberRows_ + numberColumns_);
419  lowerSlack_ = ClpCopyOfArray(rhs.lowerSlack_, numberRows_ + numberColumns_);
420  diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_ + numberColumns_);
421  deltaX_ = ClpCopyOfArray(rhs.deltaX_, numberRows_ + numberColumns_);
422  deltaZ_ = ClpCopyOfArray(rhs.deltaZ_, numberRows_ + numberColumns_);
423  deltaW_ = ClpCopyOfArray(rhs.deltaW_, numberRows_ + numberColumns_);
424  deltaSU_ = ClpCopyOfArray(rhs.deltaSU_, numberRows_ + numberColumns_);
425  deltaSL_ = ClpCopyOfArray(rhs.deltaSL_, numberRows_ + numberColumns_);
426  primalR_ = ClpCopyOfArray(rhs.primalR_, numberRows_ + numberColumns_);
427  dualR_ = ClpCopyOfArray(rhs.dualR_, numberRows_ + numberColumns_);
428  rhsB_ = ClpCopyOfArray(rhs.rhsB_, numberRows_);
429  rhsU_ = ClpCopyOfArray(rhs.rhsU_, numberRows_ + numberColumns_);
430  rhsL_ = ClpCopyOfArray(rhs.rhsL_, numberRows_ + numberColumns_);
431  rhsZ_ = ClpCopyOfArray(rhs.rhsZ_, numberRows_ + numberColumns_);
432  rhsW_ = ClpCopyOfArray(rhs.rhsW_, numberRows_ + numberColumns_);
433  rhsC_ = ClpCopyOfArray(rhs.rhsC_, numberRows_ + numberColumns_);
434  solution_ = ClpCopyOfArray(rhs.solution_, numberRows_ + numberColumns_);
435  workArray_ = ClpCopyOfArray(rhs.workArray_, numberRows_ + numberColumns_);
436  zVec_ = ClpCopyOfArray(rhs.zVec_, numberRows_ + numberColumns_);
437  wVec_ = ClpCopyOfArray(rhs.wVec_, numberRows_ + numberColumns_);
438  cholesky_ = rhs.cholesky_->clone();
439  numberComplementarityPairs_ = rhs.numberComplementarityPairs_;
440  numberComplementarityItems_ = rhs.numberComplementarityItems_;
441  maximumBarrierIterations_ = rhs.maximumBarrierIterations_;
442  gonePrimalFeasible_ = rhs.gonePrimalFeasible_;
443  goneDualFeasible_ = rhs.goneDualFeasible_;
444  algorithm_ = rhs.algorithm_;
445}
446
447void ClpInterior::gutsOfDelete()
448{
449  delete[] lower_;
450  lower_ = NULL;
451  rowLowerWork_ = NULL;
452  columnLowerWork_ = NULL;
453  delete[] upper_;
454  upper_ = NULL;
455  rowUpperWork_ = NULL;
456  columnUpperWork_ = NULL;
457  delete[] cost_;
458  cost_ = NULL;
459  delete[] rhs_;
460  rhs_ = NULL;
461  delete[] x_;
462  x_ = NULL;
463  delete[] y_;
464  y_ = NULL;
465  delete[] dj_;
466  dj_ = NULL;
467  delete lsqrObject_;
468  lsqrObject_ = NULL;
469  //delete pdcoStuff_;  // FIXME
470  pdcoStuff_ = NULL;
471  delete[] errorRegion_;
472  errorRegion_ = NULL;
473  delete[] rhsFixRegion_;
474  rhsFixRegion_ = NULL;
475  delete[] deltaY_;
476  deltaY_ = NULL;
477  delete[] upperSlack_;
478  upperSlack_ = NULL;
479  delete[] lowerSlack_;
480  lowerSlack_ = NULL;
481  delete[] diagonal_;
482  diagonal_ = NULL;
483  delete[] deltaX_;
484  deltaX_ = NULL;
485  delete[] deltaZ_;
486  deltaZ_ = NULL;
487  delete[] deltaW_;
488  deltaW_ = NULL;
489  delete[] deltaSU_;
490  deltaSU_ = NULL;
491  delete[] deltaSL_;
492  deltaSL_ = NULL;
493  delete[] primalR_;
494  primalR_ = NULL;
495  delete[] dualR_;
496  dualR_ = NULL;
497  delete[] rhsB_;
498  rhsB_ = NULL;
499  delete[] rhsU_;
500  rhsU_ = NULL;
501  delete[] rhsL_;
502  rhsL_ = NULL;
503  delete[] rhsZ_;
504  rhsZ_ = NULL;
505  delete[] rhsW_;
506  rhsW_ = NULL;
507  delete[] rhsC_;
508  rhsC_ = NULL;
509  delete[] solution_;
510  solution_ = NULL;
511  delete[] workArray_;
512  workArray_ = NULL;
513  delete[] zVec_;
514  zVec_ = NULL;
515  delete[] wVec_;
516  wVec_ = NULL;
517  delete cholesky_;
518}
519bool ClpInterior::createWorkingData()
520{
521  bool goodMatrix = true;
522  //check matrix
523  if (!matrix_->allElementsInRange(this, 1.0e-12, 1.0e20)) {
524    problemStatus_ = 4;
525    goodMatrix = false;
526  }
527  int nTotal = numberRows_ + numberColumns_;
528  delete[] solution_;
529  solution_ = new CoinWorkDouble[nTotal];
530  CoinMemcpyN(columnActivity_, numberColumns_, solution_);
531  CoinMemcpyN(rowActivity_, numberRows_, solution_ + numberColumns_);
532  delete[] cost_;
533  cost_ = new CoinWorkDouble[nTotal];
534  int i;
535  CoinWorkDouble direction = optimizationDirection_ * objectiveScale_;
536  // direction is actually scale out not scale in
537  if (direction)
538    direction = 1.0 / direction;
539  const double *obj = objective();
540  for (i = 0; i < numberColumns_; i++)
541    cost_[i] = direction * obj[i];
542  memset(cost_ + numberColumns_, 0, numberRows_ * sizeof(CoinWorkDouble));
543  // do scaling if needed
544  if (scalingFlag_ > 0 && !rowScale_) {
545    if (matrix_->scale(this))
546      scalingFlag_ = -scalingFlag_; // not scaled after all
547  }
548  delete[] lower_;
549  delete[] upper_;
550  lower_ = new CoinWorkDouble[nTotal];
551  upper_ = new CoinWorkDouble[nTotal];
552  rowLowerWork_ = lower_ + numberColumns_;
553  columnLowerWork_ = lower_;
554  rowUpperWork_ = upper_ + numberColumns_;
555  columnUpperWork_ = upper_;
556  CoinMemcpyN(rowLower_, numberRows_, rowLowerWork_);
557  CoinMemcpyN(rowUpper_, numberRows_, rowUpperWork_);
558  CoinMemcpyN(columnLower_, numberColumns_, columnLowerWork_);
559  CoinMemcpyN(columnUpper_, numberColumns_, columnUpperWork_);
560  // clean up any mismatches on infinity
561  for (i = 0; i < numberColumns_; i++) {
562    if (columnLowerWork_[i] < -1.0e30)
563      columnLowerWork_[i] = -COIN_DBL_MAX;
564    if (columnUpperWork_[i] > 1.0e30)
565      columnUpperWork_[i] = COIN_DBL_MAX;
566  }
567  // clean up any mismatches on infinity
568  for (i = 0; i < numberRows_; i++) {
569    if (rowLowerWork_[i] < -1.0e30)
570      rowLowerWork_[i] = -COIN_DBL_MAX;
571    if (rowUpperWork_[i] > 1.0e30)
572      rowUpperWork_[i] = COIN_DBL_MAX;
573  }
574  // check rim of problem okay
575  if (!sanityCheck())
576    goodMatrix = false;
577  if (rowScale_) {
578    for (i = 0; i < numberColumns_; i++) {
579      CoinWorkDouble multiplier = rhsScale_ / columnScale_[i];
580      cost_[i] *= columnScale_[i];
581      if (columnLowerWork_[i] > -1.0e50)
582        columnLowerWork_[i] *= multiplier;
583      if (columnUpperWork_[i] < 1.0e50)
584        columnUpperWork_[i] *= multiplier;
585    }
586    for (i = 0; i < numberRows_; i++) {
587      CoinWorkDouble multiplier = rhsScale_ * rowScale_[i];
588      if (rowLowerWork_[i] > -1.0e50)
589        rowLowerWork_[i] *= multiplier;
590      if (rowUpperWork_[i] < 1.0e50)
591        rowUpperWork_[i] *= multiplier;
592    }
593  } else if (rhsScale_ != 1.0) {
594    for (i = 0; i < numberColumns_ + numberRows_; i++) {
595      if (lower_[i] > -1.0e50)
596        lower_[i] *= rhsScale_;
597      if (upper_[i] < 1.0e50)
598        upper_[i] *= rhsScale_;
599    }
600  }
601  assert(!errorRegion_);
602  errorRegion_ = new CoinWorkDouble[numberRows_];
603  assert(!rhsFixRegion_);
604  rhsFixRegion_ = new CoinWorkDouble[numberRows_];
605  assert(!deltaY_);
606  deltaY_ = new CoinWorkDouble[numberRows_];
607  CoinZeroN(deltaY_, numberRows_);
608  assert(!upperSlack_);
609  upperSlack_ = new CoinWorkDouble[nTotal];
610  assert(!lowerSlack_);
611  lowerSlack_ = new CoinWorkDouble[nTotal];
612  assert(!diagonal_);
613  diagonal_ = new CoinWorkDouble[nTotal];
614  assert(!deltaX_);
615  deltaX_ = new CoinWorkDouble[nTotal];
616  CoinZeroN(deltaX_, nTotal);
617  assert(!deltaZ_);
618  deltaZ_ = new CoinWorkDouble[nTotal];
619  CoinZeroN(deltaZ_, nTotal);
620  assert(!deltaW_);
621  deltaW_ = new CoinWorkDouble[nTotal];
622  CoinZeroN(deltaW_, nTotal);
623  assert(!deltaSU_);
624  deltaSU_ = new CoinWorkDouble[nTotal];
625  CoinZeroN(deltaSU_, nTotal);
626  assert(!deltaSL_);
627  deltaSL_ = new CoinWorkDouble[nTotal];
628  CoinZeroN(deltaSL_, nTotal);
629  assert(!primalR_);
630  assert(!dualR_);
631  // create arrays if we are doing KKT
632  if (cholesky_->type() >= 20) {
633    primalR_ = new CoinWorkDouble[nTotal];
634    CoinZeroN(primalR_, nTotal);
635    dualR_ = new CoinWorkDouble[numberRows_];
636    CoinZeroN(dualR_, numberRows_);
637  }
638  assert(!rhsB_);
639  rhsB_ = new CoinWorkDouble[numberRows_];
640  CoinZeroN(rhsB_, numberRows_);
641  assert(!rhsU_);
642  rhsU_ = new CoinWorkDouble[nTotal];
643  CoinZeroN(rhsU_, nTotal);
644  assert(!rhsL_);
645  rhsL_ = new CoinWorkDouble[nTotal];
646  CoinZeroN(rhsL_, nTotal);
647  assert(!rhsZ_);
648  rhsZ_ = new CoinWorkDouble[nTotal];
649  CoinZeroN(rhsZ_, nTotal);
650  assert(!rhsW_);
651  rhsW_ = new CoinWorkDouble[nTotal];
652  CoinZeroN(rhsW_, nTotal);
653  assert(!rhsC_);
654  rhsC_ = new CoinWorkDouble[nTotal];
655  CoinZeroN(rhsC_, nTotal);
656  assert(!workArray_);
657  workArray_ = new CoinWorkDouble[nTotal];
658  CoinZeroN(workArray_, nTotal);
659  assert(!zVec_);
660  zVec_ = new CoinWorkDouble[nTotal];
661  CoinZeroN(zVec_, nTotal);
662  assert(!wVec_);
663  wVec_ = new CoinWorkDouble[nTotal];
664  CoinZeroN(wVec_, nTotal);
665  assert(!dj_);
666  dj_ = new CoinWorkDouble[nTotal];
667  if (!status_)
668    status_ = new unsigned char[numberRows_ + numberColumns_];
669  memset(status_, 0, numberRows_ + numberColumns_);
670  return goodMatrix;
671}
672void ClpInterior::deleteWorkingData()
673{
674  int i;
675  if (optimizationDirection_ != 1.0 || objectiveScale_ != 1.0) {
676    CoinWorkDouble scaleC = optimizationDirection_ / objectiveScale_;
677    // and modify all dual signs
678    for (i = 0; i < numberColumns_; i++)
679      reducedCost_[i] = scaleC * dj_[i];
680    for (i = 0; i < numberRows_; i++)
681      dual_[i] *= scaleC;
682  }
683  if (rowScale_) {
684    CoinWorkDouble scaleR = 1.0 / rhsScale_;
685    for (i = 0; i < numberColumns_; i++) {
686      CoinWorkDouble scaleFactor = columnScale_[i];
687      CoinWorkDouble valueScaled = columnActivity_[i];
688      columnActivity_[i] = valueScaled * scaleFactor * scaleR;
689      CoinWorkDouble valueScaledDual = reducedCost_[i];
690      reducedCost_[i] = valueScaledDual / scaleFactor;
691    }
692    for (i = 0; i < numberRows_; i++) {
693      CoinWorkDouble scaleFactor = rowScale_[i];
694      CoinWorkDouble valueScaled = rowActivity_[i];
695      rowActivity_[i] = (valueScaled * scaleR) / scaleFactor;
696      CoinWorkDouble valueScaledDual = dual_[i];
697      dual_[i] = valueScaledDual * scaleFactor;
698    }
699  } else if (rhsScale_ != 1.0) {
700    CoinWorkDouble scaleR = 1.0 / rhsScale_;
701    for (i = 0; i < numberColumns_; i++) {
702      CoinWorkDouble valueScaled = columnActivity_[i];
703      columnActivity_[i] = valueScaled * scaleR;
704    }
705    for (i = 0; i < numberRows_; i++) {
706      CoinWorkDouble valueScaled = rowActivity_[i];
707      rowActivity_[i] = valueScaled * scaleR;
708    }
709  }
710  delete[] cost_;
711  cost_ = NULL;
712  delete[] solution_;
713  solution_ = NULL;
714  delete[] lower_;
715  lower_ = NULL;
716  delete[] upper_;
717  upper_ = NULL;
718  delete[] errorRegion_;
719  errorRegion_ = NULL;
720  delete[] rhsFixRegion_;
721  rhsFixRegion_ = NULL;
722  delete[] deltaY_;
723  deltaY_ = NULL;
724  delete[] upperSlack_;
725  upperSlack_ = NULL;
726  delete[] lowerSlack_;
727  lowerSlack_ = NULL;
728  delete[] diagonal_;
729  diagonal_ = NULL;
730  delete[] deltaX_;
731  deltaX_ = NULL;
732  delete[] workArray_;
733  workArray_ = NULL;
734  delete[] zVec_;
735  zVec_ = NULL;
736  delete[] wVec_;
737  wVec_ = NULL;
738  delete[] dj_;
739  dj_ = NULL;
740}
741// Sanity check on input data - returns true if okay
742bool ClpInterior::sanityCheck()
743{
744  // bad if empty
745  if (!numberColumns_ || ((!numberRows_ || !matrix_->getNumElements()) && objective_->type() < 2)) {
746    problemStatus_ = emptyProblem();
747    return false;
748  }
749  int numberBad;
750  CoinWorkDouble largestBound, smallestBound, minimumGap;
751  CoinWorkDouble smallestObj, largestObj;
752  int firstBad;
753  int modifiedBounds = 0;
754  int i;
755  numberBad = 0;
756  firstBad = -1;
757  minimumGap = 1.0e100;
758  smallestBound = 1.0e100;
759  largestBound = 0.0;
760  smallestObj = 1.0e100;
761  largestObj = 0.0;
762  // If bounds are too close - fix
763  CoinWorkDouble fixTolerance = 1.1 * primalTolerance();
764  for (i = numberColumns_; i < numberColumns_ + numberRows_; i++) {
765    CoinWorkDouble value;
766    value = CoinAbs(cost_[i]);
767    if (value > 1.0e50) {
768      numberBad++;
769      if (firstBad < 0)
770        firstBad = i;
771    } else if (value) {
772      if (value > largestObj)
773        largestObj = value;
774      if (value < smallestObj)
775        smallestObj = value;
776    }
777    value = upper_[i] - lower_[i];
778    if (value < -primalTolerance()) {
779      numberBad++;
780      if (firstBad < 0)
781        firstBad = i;
782    } else if (value <= fixTolerance) {
783      if (value) {
784        // modify
785        upper_[i] = lower_[i];
786        modifiedBounds++;
787      }
788    } else {
789      if (value < minimumGap)
790        minimumGap = value;
791    }
792    if (lower_[i] > -1.0e100 && lower_[i]) {
793      value = CoinAbs(lower_[i]);
794      if (value > largestBound)
795        largestBound = value;
796      if (value < smallestBound)
797        smallestBound = value;
798    }
799    if (upper_[i] < 1.0e100 && upper_[i]) {
800      value = CoinAbs(upper_[i]);
801      if (value > largestBound)
802        largestBound = value;
803      if (value < smallestBound)
804        smallestBound = value;
805    }
806  }
807  if (largestBound)
808    handler_->message(CLP_RIMSTATISTICS3, messages_)
809      << static_cast< double >(smallestBound)
810      << static_cast< double >(largestBound)
811      << static_cast< double >(minimumGap)
812      << CoinMessageEol;
813  minimumGap = 1.0e100;
814  smallestBound = 1.0e100;
815  largestBound = 0.0;
816  for (i = 0; i < numberColumns_; i++) {
817    CoinWorkDouble value;
818    value = CoinAbs(cost_[i]);
819    if (value > 1.0e50) {
820      numberBad++;
821      if (firstBad < 0)
822        firstBad = i;
823    } else if (value) {
824      if (value > largestObj)
825        largestObj = value;
826      if (value < smallestObj)
827        smallestObj = value;
828    }
829    value = upper_[i] - lower_[i];
830    if (value < -primalTolerance()) {
831      numberBad++;
832      if (firstBad < 0)
833        firstBad = i;
834    } else if (value <= fixTolerance) {
835      if (value) {
836        // modify
837        upper_[i] = lower_[i];
838        modifiedBounds++;
839      }
840    } else {
841      if (value < minimumGap)
842        minimumGap = value;
843    }
844    if (lower_[i] > -1.0e100 && lower_[i]) {
845      value = CoinAbs(lower_[i]);
846      if (value > largestBound)
847        largestBound = value;
848      if (value < smallestBound)
849        smallestBound = value;
850    }
851    if (upper_[i] < 1.0e100 && upper_[i]) {
852      value = CoinAbs(upper_[i]);
853      if (value > largestBound)
854        largestBound = value;
855      if (value < smallestBound)
856        smallestBound = value;
857    }
858  }
859  char rowcol[] = { 'R', 'C' };
860  if (numberBad) {
861    handler_->message(CLP_BAD_BOUNDS, messages_)
862      << numberBad
863      << rowcol[isColumn(firstBad)] << sequenceWithin(firstBad)
864      << CoinMessageEol;
865    problemStatus_ = 4;
866    return false;
867  }
868  if (modifiedBounds)
869    handler_->message(CLP_MODIFIEDBOUNDS, messages_)
870      << modifiedBounds
871      << CoinMessageEol;
872  handler_->message(CLP_RIMSTATISTICS1, messages_)
873    << static_cast< double >(smallestObj)
874    << static_cast< double >(largestObj)
875    << CoinMessageEol;
876  if (largestBound)
877    handler_->message(CLP_RIMSTATISTICS2, messages_)
878      << static_cast< double >(smallestBound)
879      << static_cast< double >(largestBound)
880      << static_cast< double >(minimumGap)
881      << CoinMessageEol;
882  return true;
883}
884/* Loads a problem (the constraints on the
885   rows are given by lower and upper bounds). If a pointer is 0 then the
886   following values are the default:
887   <ul>
888   <li> <code>colub</code>: all columns have upper bound infinity
889   <li> <code>collb</code>: all columns have lower bound 0
890   <li> <code>rowub</code>: all rows have upper bound infinity
891   <li> <code>rowlb</code>: all rows have lower bound -infinity
892   <li> <code>obj</code>: all variables have 0 objective coefficient
893   </ul>
894*/
895void ClpInterior::loadProblem(const ClpMatrixBase &matrix,
896  const double *collb, const double *colub,
897  const double *obj,
898  const double *rowlb, const double *rowub,
899  const double *rowObjective)
900{
901  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
902    rowObjective);
903}
904void ClpInterior::loadProblem(const CoinPackedMatrix &matrix,
905  const double *collb, const double *colub,
906  const double *obj,
907  const double *rowlb, const double *rowub,
908  const double *rowObjective)
909{
910  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
911    rowObjective);
912}
913
914/* Just like the other loadProblem() method except that the matrix is
915   given in a standard column major ordered format (without gaps). */
916void ClpInterior::loadProblem(const int numcols, const int numrows,
917  const CoinBigIndex *start, const int *index,
918  const double *value,
919  const double *collb, const double *colub,
920  const double *obj,
921  const double *rowlb, const double *rowub,
922  const double *rowObjective)
923{
924  ClpModel::loadProblem(numcols, numrows, start, index, value,
925    collb, colub, obj, rowlb, rowub,
926    rowObjective);
927}
928void ClpInterior::loadProblem(const int numcols, const int numrows,
929  const CoinBigIndex *start, const int *index,
930  const double *value, const int *length,
931  const double *collb, const double *colub,
932  const double *obj,
933  const double *rowlb, const double *rowub,
934  const double *rowObjective)
935{
936  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
937    collb, colub, obj, rowlb, rowub,
938    rowObjective);
939}
940// Read an mps file from the given filename
941int ClpInterior::readMps(const char *filename,
942  bool keepNames,
943  bool ignoreErrors)
944{
945  int status = ClpModel::readMps(filename, keepNames, ignoreErrors);
946  return status;
947}
948#include "ClpPdco.hpp"
949/* Pdco algorithm - see ClpPdco.hpp for method */
950int ClpInterior::pdco()
951{
952  return ((ClpPdco *)this)->pdco();
953}
954// ** Temporary version
955int ClpInterior::pdco(ClpPdcoBase *stuff, Options &options, Info &info, Outfo &outfo)
956{
957  return ((ClpPdco *)this)->pdco(stuff, options, info, outfo);
958}
959#include "ClpPredictorCorrector.hpp"
960// Primal-Dual Predictor-Corrector barrier
961int ClpInterior::primalDual()
962{
963  return (static_cast< ClpPredictorCorrector * >(this))->solve();
964}
965
966void ClpInterior::checkSolution()
967{
968  int iRow, iColumn;
969  CoinWorkDouble *reducedCost = reinterpret_cast< CoinWorkDouble * >(reducedCost_);
970  CoinWorkDouble *dual = reinterpret_cast< CoinWorkDouble * >(dual_);
971  CoinMemcpyN(cost_, numberColumns_, reducedCost);
972  matrix_->transposeTimes(-1.0, dual, reducedCost);
973  // Now modify reduced costs for quadratic
974  CoinWorkDouble quadraticOffset = quadraticDjs(reducedCost,
975    solution_, scaleFactor_);
976
977  objectiveValue_ = 0.0;
978  // now look at solution
979  sumPrimalInfeasibilities_ = 0.0;
980  sumDualInfeasibilities_ = 0.0;
981  CoinWorkDouble dualTolerance = 10.0 * dblParam_[ClpDualTolerance];
982  CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance];
983  CoinWorkDouble primalTolerance2 = 10.0 * dblParam_[ClpPrimalTolerance];
984  worstComplementarity_ = 0.0;
985  complementarityGap_ = 0.0;
986
987  // Done scaled - use permanent regions for output
988  // but internal for bounds
989  const CoinWorkDouble *lower = lower_ + numberColumns_;
990  const CoinWorkDouble *upper = upper_ + numberColumns_;
991  for (iRow = 0; iRow < numberRows_; iRow++) {
992    CoinWorkDouble infeasibility = 0.0;
993    CoinWorkDouble distanceUp = CoinMin(upper[iRow] - rowActivity_[iRow], static_cast< CoinWorkDouble >(1.0e10));
994    CoinWorkDouble distanceDown = CoinMin(rowActivity_[iRow] - lower[iRow], static_cast< CoinWorkDouble >(1.0e10));
995    if (distanceUp > primalTolerance2) {
996      CoinWorkDouble value = dual[iRow];
997      // should not be negative
998      if (value < -dualTolerance) {
999        sumDualInfeasibilities_ += -dualTolerance - value;
1000        value = -value * distanceUp;
1001        if (value > worstComplementarity_)
1002          worstComplementarity_ = value;
1003        complementarityGap_ += value;
1004      }
1005    }
1006    if (distanceDown > primalTolerance2) {
1007      CoinWorkDouble value = dual[iRow];
1008      // should not be positive
1009      if (value > dualTolerance) {
1010        sumDualInfeasibilities_ += value - dualTolerance;
1011        value = value * distanceDown;
1012        if (value > worstComplementarity_)
1013          worstComplementarity_ = value;
1014        complementarityGap_ += value;
1015      }
1016    }
1017    if (rowActivity_[iRow] > upper[iRow]) {
1018      infeasibility = rowActivity_[iRow] - upper[iRow];
1019    } else if (rowActivity_[iRow] < lower[iRow]) {
1020      infeasibility = lower[iRow] - rowActivity_[iRow];
1021    }
1022    if (infeasibility > primalTolerance) {
1023      sumPrimalInfeasibilities_ += infeasibility - primalTolerance;
1024    }
1025  }
1026  lower = lower_;
1027  upper = upper_;
1028  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1029    CoinWorkDouble infeasibility = 0.0;
1030    objectiveValue_ += cost_[iColumn] * columnActivity_[iColumn];
1031    CoinWorkDouble distanceUp = CoinMin(upper[iColumn] - columnActivity_[iColumn], static_cast< CoinWorkDouble >(1.0e10));
1032    CoinWorkDouble distanceDown = CoinMin(columnActivity_[iColumn] - lower[iColumn], static_cast< CoinWorkDouble >(1.0e10));
1033    if (distanceUp > primalTolerance2) {
1034      CoinWorkDouble value = reducedCost[iColumn];
1035      // should not be negative
1036      if (value < -dualTolerance) {
1037        sumDualInfeasibilities_ += -dualTolerance - value;
1038        value = -value * distanceUp;
1039        if (value > worstComplementarity_)
1040          worstComplementarity_ = value;
1041        complementarityGap_ += value;
1042      }
1043    }
1044    if (distanceDown > primalTolerance2) {
1045      CoinWorkDouble value = reducedCost[iColumn];
1046      // should not be positive
1047      if (value > dualTolerance) {
1048        sumDualInfeasibilities_ += value - dualTolerance;
1049        value = value * distanceDown;
1050        if (value > worstComplementarity_)
1051          worstComplementarity_ = value;
1052        complementarityGap_ += value;
1053      }
1054    }
1055    if (columnActivity_[iColumn] > upper[iColumn]) {
1056      infeasibility = columnActivity_[iColumn] - upper[iColumn];
1057    } else if (columnActivity_[iColumn] < lower[iColumn]) {
1058      infeasibility = lower[iColumn] - columnActivity_[iColumn];
1059    }
1060    if (infeasibility > primalTolerance) {
1061      sumPrimalInfeasibilities_ += infeasibility - primalTolerance;
1062    }
1063  }
1064#if COIN_LONG_WORK
1065  // ok as packs down
1066  CoinMemcpyN(reducedCost, numberColumns_, reducedCost_);
1067#endif
1068  // add in offset
1069  objectiveValue_ += 0.5 * quadraticOffset;
1070}
1071// Set cholesky (and delete present one)
1072void ClpInterior::setCholesky(ClpCholeskyBase *cholesky)
1073{
1074  delete cholesky_;
1075  cholesky_ = cholesky;
1076}
1077/* Borrow model.  This is so we dont have to copy large amounts
1078   of data around.  It assumes a derived class wants to overwrite
1079   an empty model with a real one - while it does an algorithm.
1080   This is same as ClpModel one. */
1081void ClpInterior::borrowModel(ClpModel &otherModel)
1082{
1083  ClpModel::borrowModel(otherModel);
1084}
1085/* Return model - updates any scalars */
1086void ClpInterior::returnModel(ClpModel &otherModel)
1087{
1088  ClpModel::returnModel(otherModel);
1089}
1090// Return number fixed to see if worth presolving
1091int ClpInterior::numberFixed() const
1092{
1093  int i;
1094  int nFixed = 0;
1095  for (i = 0; i < numberColumns_; i++) {
1096    if (columnUpper_[i] < 1.0e20 || columnLower_[i] > -1.0e20) {
1097      if (columnUpper_[i] > columnLower_[i]) {
1098        if (fixedOrFree(i))
1099          nFixed++;
1100      }
1101    }
1102  }
1103  for (i = 0; i < numberRows_; i++) {
1104    if (rowUpper_[i] < 1.0e20 || rowLower_[i] > -1.0e20) {
1105      if (rowUpper_[i] > rowLower_[i]) {
1106        if (fixedOrFree(i + numberColumns_))
1107          nFixed++;
1108      }
1109    }
1110  }
1111  return nFixed;
1112}
1113// fix variables interior says should be
1114void ClpInterior::fixFixed(bool reallyFix)
1115{
1116  // Arrays for change in columns and rhs
1117  CoinWorkDouble *columnChange = new CoinWorkDouble[numberColumns_];
1118  CoinWorkDouble *rowChange = new CoinWorkDouble[numberRows_];
1119  CoinZeroN(columnChange, numberColumns_);
1120  CoinZeroN(rowChange, numberRows_);
1121  matrix_->times(1.0, columnChange, rowChange);
1122  int i;
1123  CoinWorkDouble tolerance = primalTolerance();
1124  for (i = 0; i < numberColumns_; i++) {
1125    if (columnUpper_[i] < 1.0e20 || columnLower_[i] > -1.0e20) {
1126      if (columnUpper_[i] > columnLower_[i]) {
1127        if (fixedOrFree(i)) {
1128          if (columnActivity_[i] - columnLower_[i] < columnUpper_[i] - columnActivity_[i]) {
1129            CoinWorkDouble change = columnLower_[i] - columnActivity_[i];
1130            if (CoinAbs(change) < tolerance) {
1131              if (reallyFix)
1132                columnUpper_[i] = columnLower_[i];
1133              columnChange[i] = change;
1134              columnActivity_[i] = columnLower_[i];
1135            }
1136          } else {
1137            CoinWorkDouble change = columnUpper_[i] - columnActivity_[i];
1138            if (CoinAbs(change) < tolerance) {
1139              if (reallyFix)
1140                columnLower_[i] = columnUpper_[i];
1141              columnChange[i] = change;
1142              columnActivity_[i] = columnUpper_[i];
1143            }
1144          }
1145        }
1146      }
1147    }
1148  }
1149  CoinZeroN(rowChange, numberRows_);
1150  matrix_->times(1.0, columnChange, rowChange);
1151  // If makes mess of things then don't do
1152  CoinWorkDouble newSum = 0.0;
1153  for (i = 0; i < numberRows_; i++) {
1154    CoinWorkDouble value = rowActivity_[i] + rowChange[i];
1155    if (value > rowUpper_[i] + tolerance)
1156      newSum += value - rowUpper_[i] - tolerance;
1157    else if (value < rowLower_[i] - tolerance)
1158      newSum -= value - rowLower_[i] + tolerance;
1159  }
1160  if (newSum > 1.0e-5 + 1.5 * sumPrimalInfeasibilities_) {
1161    // put back and skip changes
1162    for (i = 0; i < numberColumns_; i++)
1163      columnActivity_[i] -= columnChange[i];
1164  } else {
1165    CoinZeroN(rowActivity_, numberRows_);
1166    matrix_->times(1.0, columnActivity_, rowActivity_);
1167    if (reallyFix) {
1168      for (i = 0; i < numberRows_; i++) {
1169        if (rowUpper_[i] < 1.0e20 || rowLower_[i] > -1.0e20) {
1170          if (rowUpper_[i] > rowLower_[i]) {
1171            if (fixedOrFree(i + numberColumns_)) {
1172              if (rowActivity_[i] - rowLower_[i] < rowUpper_[i] - rowActivity_[i]) {
1173                CoinWorkDouble change = rowLower_[i] - rowActivity_[i];
1174                if (CoinAbs(change) < tolerance) {
1175                  if (reallyFix)
1176                    rowUpper_[i] = rowLower_[i];
1177                  rowActivity_[i] = rowLower_[i];
1178                }
1179              } else {
1180                CoinWorkDouble change = rowLower_[i] - rowActivity_[i];
1181                if (CoinAbs(change) < tolerance) {
1182                  if (reallyFix)
1183                    rowLower_[i] = rowUpper_[i];
1184                  rowActivity_[i] = rowUpper_[i];
1185                }
1186              }
1187            }
1188          }
1189        }
1190      }
1191    }
1192  }
1193  delete[] rowChange;
1194  delete[] columnChange;
1195}
1196/* Modifies djs to allow for quadratic.
1197   returns quadratic offset */
1198CoinWorkDouble
1199ClpInterior::quadraticDjs(CoinWorkDouble *djRegion, const CoinWorkDouble *solution,
1200  CoinWorkDouble scaleFactor)
1201{
1202  CoinWorkDouble quadraticOffset = 0.0;
1203#ifndef NO_RTTI
1204  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_));
1205#else
1206  ClpQuadraticObjective *quadraticObj = NULL;
1207  if (objective_->type() == 2)
1208    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
1209#endif
1210  if (quadraticObj) {
1211    CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
1212    const int *columnQuadratic = quadratic->getIndices();
1213    const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
1214    const int *columnQuadraticLength = quadratic->getVectorLengths();
1215    double *quadraticElement = quadratic->getMutableElements();
1216    int numberColumns = quadratic->getNumCols();
1217    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1218      CoinWorkDouble value = 0.0;
1219      for (CoinBigIndex j = columnQuadraticStart[iColumn];
1220           j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1221        int jColumn = columnQuadratic[j];
1222        CoinWorkDouble valueJ = solution[jColumn];
1223        CoinWorkDouble elementValue = quadraticElement[j];
1224        //value += valueI*valueJ*elementValue;
1225        value += valueJ * elementValue;
1226        quadraticOffset += solution[iColumn] * valueJ * elementValue;
1227      }
1228      djRegion[iColumn] += scaleFactor * value;
1229    }
1230  }
1231  return quadraticOffset;
1232}
1233
1234/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1235*/
Note: See TracBrowser for help on using the repository browser.