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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 43.5 KB
Line 
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}
Note: See TracBrowser for help on using the repository browser.