source: stable/1.15/Clp/src/ClpInterior.cpp @ 1941

Last change on this file since 1941 was 1941, checked in by stefan, 6 years ago

sync with trunk rev 1940

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