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