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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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