1/* $Id: ClpFactorization.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include "ClpFactorization.hpp"
8#ifndef SLIM_CLP
9#include "ClpQuadraticObjective.hpp"
10#endif
11#include "CoinHelperFunctions.hpp"
12#include "CoinIndexedVector.hpp"
13#include "ClpSimplex.hpp"
14#include "ClpSimplexDual.hpp"
15#include "ClpMatrixBase.hpp"
16#ifndef SLIM_CLP
17#include "ClpNetworkBasis.hpp"
18#include "ClpNetworkMatrix.hpp"
19//#define CHECK_NETWORK
20#ifdef CHECK_NETWORK
21const static bool doCheck = true;
22#else
23const static bool doCheck = false;
24#endif
25#endif
26//#define CLP_FACTORIZATION_INSTRUMENT
27#ifdef CLP_FACTORIZATION_INSTRUMENT
28#include "CoinTime.hpp"
29double factorization_instrument(int type)
30{
31 static int times[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
32 static double startTime = 0.0;
33 static double totalTimes [10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
34 if (type < 0) {
35 assert (!startTime);
36 startTime = CoinCpuTime();
37 return 0.0;
38 } else if (type > 0) {
39 times[type]++;
40 double difference = CoinCpuTime() - startTime;
41 totalTimes[type] += difference;
42 startTime = 0.0;
43 return difference;
44 } else {
45 // report
46 const char * types[10] = {
47 "", "fac=rhs_etc", "factorize", "replace", "update_FT",
48 "update", "update_transpose", "gosparse", "getWeights!", "update2_FT"
49 };
50 double total = 0.0;
51 for (int i = 1; i < 10; i++) {
52 if (times[i]) {
53 printf("%s was called %d times taking %g seconds\n",
54 types[i], times[i], totalTimes[i]);
55 total += totalTimes[i];
56 times[i] = 0;
57 totalTimes[i] = 0.0;
58 }
59 }
60 return total;
61 }
62}
63#endif
64//#############################################################################
65// Constructors / Destructor / Assignment
66//#############################################################################
67#ifndef CLP_MULTIPLE_FACTORIZATIONS
68
69//-------------------------------------------------------------------
70// Default Constructor
71//-------------------------------------------------------------------
72ClpFactorization::ClpFactorization () :
73 CoinFactorization()
74{
75#ifndef SLIM_CLP
76 networkBasis_ = NULL;
77#endif
78}
79
80//-------------------------------------------------------------------
81// Copy constructor
82//-------------------------------------------------------------------
83ClpFactorization::ClpFactorization (const ClpFactorization & rhs,
84 int dummyDenseIfSmaller) :
85 CoinFactorization(rhs)
86{
87#ifndef SLIM_CLP
88 if (rhs.networkBasis_)
89 networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
90 else
91 networkBasis_ = NULL;
92#endif
93}
94
95ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
96 CoinFactorization(rhs)
97{
98#ifndef SLIM_CLP
99 networkBasis_ = NULL;
100#endif
101}
102
103//-------------------------------------------------------------------
104// Destructor
105//-------------------------------------------------------------------
106ClpFactorization::~ClpFactorization ()
107{
108#ifndef SLIM_CLP
109 delete networkBasis_;
110#endif
111}
112
113//----------------------------------------------------------------
114// Assignment operator
115//-------------------------------------------------------------------
116ClpFactorization &
117ClpFactorization::operator=(const ClpFactorization& rhs)
118{
119#ifdef CLP_FACTORIZATION_INSTRUMENT
120 factorization_instrument(-1);
121#endif
122 if (this != &rhs) {
123 CoinFactorization::operator=(rhs);
124#ifndef SLIM_CLP
125 delete networkBasis_;
126 if (rhs.networkBasis_)
127 networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
128 else
129 networkBasis_ = NULL;
130#endif
131 }
132#ifdef CLP_FACTORIZATION_INSTRUMENT
133 factorization_instrument(1);
134#endif
135 return *this;
136}
137#if 0
138static unsigned int saveList[10000];
139int numberSave = -1;
140inline bool isDense(int i)
141{
142 return ((saveList[i>>5] >> (i & 31)) & 1) != 0;
143}
144inline void setDense(int i)
145{
146 unsigned int & value = saveList[i>>5];
147 int bit = i & 31;
148 value |= (1 << bit);
149}
150#endif
151int
152ClpFactorization::factorize ( ClpSimplex * model,
153 int solveType, bool valuesPass)
154{
155#ifdef CLP_FACTORIZATION_INSTRUMENT
156 factorization_instrument(-1);
157#endif
158 ClpMatrixBase * matrix = model->clpMatrix();
159 int numberRows = model->numberRows();
160 int numberColumns = model->numberColumns();
161 if (!numberRows)
162 return 0;
163 // If too many compressions increase area
164 if (numberPivots_ > 1 && numberCompressions_ * 10 > numberPivots_ + 10) {
165 areaFactor_ *= 1.1;
166 }
167 //int numberPivots=numberPivots_;
168#if 0
169 if (model->algorithm() > 0)
170 numberSave = -1;
171#endif
172#ifndef SLIM_CLP
173 if (!networkBasis_ || doCheck) {
174#endif
175 status_ = -99;
176 int * pivotVariable = model->pivotVariable();
177 //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
178 while (status_ < -98) {
179
180 int i;
181 int numberBasic = 0;
182 int numberRowBasic;
183 // Move pivot variables across if they look good
184 int * pivotTemp = model->rowArray(0)->getIndices();
185 assert (!model->rowArray(0)->getNumElements());
186 if (!matrix->rhsOffset(model)) {
187#if 0
188 if (numberSave > 0) {
189 int nStill = 0;
190 int nAtBound = 0;
191 int nZeroDual = 0;
192 CoinIndexedVector * array = model->rowArray(3);
193 CoinIndexedVector * objArray = model->columnArray(1);
194 array->clear();
195 objArray->clear();
196 double * cost = model->costRegion();
197 double tolerance = model->primalTolerance();
198 double offset = 0.0;
199 for (i = 0; i < numberRows; i++) {
200 int iPivot = pivotVariable[i];
201 if (iPivot < numberColumns && isDense(iPivot)) {
202 if (model->getColumnStatus(iPivot) == ClpSimplex::basic) {
203 nStill++;
204 double value = model->solutionRegion()[iPivot];
205 double dual = model->dualRowSolution()[i];
206 double lower = model->lowerRegion()[iPivot];
207 double upper = model->upperRegion()[iPivot];
208 ClpSimplex::Status status;
209 if (fabs(value - lower) < tolerance) {
210 status = ClpSimplex::atLowerBound;
211 nAtBound++;
212 } else if (fabs(value - upper) < tolerance) {
213 nAtBound++;
214 status = ClpSimplex::atUpperBound;
215 } else if (value > lower && value < upper) {
216 status = ClpSimplex::superBasic;
217 } else {
218 status = ClpSimplex::basic;
219 }
220 if (status != ClpSimplex::basic) {
221 if (model->getRowStatus(i) != ClpSimplex::basic) {
222 model->setColumnStatus(iPivot, ClpSimplex::atLowerBound);
223 model->setRowStatus(i, ClpSimplex::basic);
224 pivotVariable[i] = i + numberColumns;
225 model->dualRowSolution()[i] = 0.0;
226 model->djRegion(0)[i] = 0.0;
227 array->add(i, dual);
228 offset += dual * model->solutionRegion(0)[i];
229 }
230 }
231 if (fabs(dual) < 1.0e-5)
232 nZeroDual++;
233 }
234 }
235 }
236 printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
237 numberSave, nStill, nAtBound, nZeroDual, offset);
238 if (array->getNumElements()) {
239 // modify costs
240 model->clpMatrix()->transposeTimes(model, 1.0, array, model->columnArray(0),
241 objArray);
242 array->clear();
243 int n = objArray->getNumElements();
244 int * indices = objArray->getIndices();
245 double * elements = objArray->denseVector();
246 for (i = 0; i < n; i++) {
247 int iColumn = indices[i];
248 cost[iColumn] -= elements[iColumn];
249 elements[iColumn] = 0.0;
250 }
251 objArray->setNumElements(0);
252 }
253 }
254#endif
255 // Seems to prefer things in order so quickest
256 // way is to go though like this
257 for (i = 0; i < numberRows; i++) {
258 if (model->getRowStatus(i) == ClpSimplex::basic)
259 pivotTemp[numberBasic++] = i;
260 }
261 numberRowBasic = numberBasic;
262 /* Put column basic variables into pivotVariable
263 This is done by ClpMatrixBase to allow override for gub
264 */
265 matrix->generalExpanded(model, 0, numberBasic);
266 } else {
267 // Long matrix - do a different way
268 bool fullSearch = false;
269 for (i = 0; i < numberRows; i++) {
270 int iPivot = pivotVariable[i];
271 if (iPivot >= numberColumns) {
272 pivotTemp[numberBasic++] = iPivot - numberColumns;
273 }
274 }
275 numberRowBasic = numberBasic;
276 for (i = 0; i < numberRows; i++) {
277 int iPivot = pivotVariable[i];
278 if (iPivot < numberColumns) {
279 if (iPivot >= 0) {
280 pivotTemp[numberBasic++] = iPivot;
281 } else {
282 // not full basis
283 fullSearch = true;
284 break;
285 }
286 }
287 }
288 if (fullSearch) {
289 // do slow way
290 numberBasic = 0;
291 for (i = 0; i < numberRows; i++) {
292 if (model->getRowStatus(i) == ClpSimplex::basic)
293 pivotTemp[numberBasic++] = i;
294 }
295 numberRowBasic = numberBasic;
296 /* Put column basic variables into pivotVariable
297 This is done by ClpMatrixBase to allow override for gub
298 */
299 matrix->generalExpanded(model, 0, numberBasic);
300 }
301 }
302 if (numberBasic > model->maximumBasic()) {
303#if 0 // ndef NDEBUG
304 printf("%d basic - should only be %d\n",
305 numberBasic, numberRows);
306#endif
307 // Take out some
308 numberBasic = numberRowBasic;
309 for (int i = 0; i < numberColumns; i++) {
310 if (model->getColumnStatus(i) == ClpSimplex::basic) {
311 if (numberBasic < numberRows)
312 numberBasic++;
313 else
314 model->setColumnStatus(i, ClpSimplex::superBasic);
315 }
316 }
317 numberBasic = numberRowBasic;
318 matrix->generalExpanded(model, 0, numberBasic);
319 }
320#ifndef SLIM_CLP
321 // see if matrix a network
322#ifndef NO_RTTI
323 ClpNetworkMatrix* networkMatrix =
324 dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
325#else
326 ClpNetworkMatrix* networkMatrix = NULL;
327 if (model->clpMatrix()->type() == 11)
328 networkMatrix =
329 static_cast< ClpNetworkMatrix*>(model->clpMatrix());
330#endif
331 // If network - still allow ordinary factorization first time for laziness
332 if (networkMatrix)
333 biasLU_ = 0; // All to U if network
334 //int saveMaximumPivots = maximumPivots();
335 delete networkBasis_;
336 networkBasis_ = NULL;
337 if (networkMatrix && !doCheck)
338 maximumPivots(1);
339#endif
340 //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
341 while (status_ == -99) {
342 // maybe for speed will be better to leave as many regions as possible
343 gutsOfDestructor();
344 gutsOfInitialize(2);
345 CoinBigIndex numberElements = numberRowBasic;
346
347 // compute how much in basis
348
349 int i;
350 // can change for gub
351 int numberColumnBasic = numberBasic - numberRowBasic;
352
353 numberElements += matrix->countBasis(model,
354 pivotTemp + numberRowBasic,
355 numberRowBasic,
356 numberColumnBasic);
357 // and recompute as network side say different
358 if (model->numberIterations())
359 numberRowBasic = numberBasic - numberColumnBasic;
360 numberElements = 3 * numberBasic + 3 * numberElements + 20000;
361#if 0
362 // If iteration not zero then may be compressed
363 getAreas ( !model->numberIterations() ? numberRows : numberBasic,
364 numberRowBasic + numberColumnBasic, numberElements,
365 2 * numberElements );
366#else
367 getAreas ( numberRows,
368 numberRowBasic + numberColumnBasic, numberElements,
369 2 * numberElements );
370#endif
371 //fill
372 // Fill in counts so we can skip part of preProcess
373 int * numberInRow = numberInRow_.array();
374 int * numberInColumn = numberInColumn_.array();
375 CoinZeroN ( numberInRow, numberRows_ + 1 );
376 CoinZeroN ( numberInColumn, maximumColumnsExtra_ + 1 );
377 double * elementU = elementU_.array();
378 int * indexRowU = indexRowU_.array();
379 CoinBigIndex * startColumnU = startColumnU_.array();
380 for (i = 0; i < numberRowBasic; i++) {
381 int iRow = pivotTemp[i];
382 indexRowU[i] = iRow;
383 startColumnU[i] = i;
384 elementU[i] = slackValue_;
385 numberInRow[iRow] = 1;
386 numberInColumn[i] = 1;
387 }
388 startColumnU[numberRowBasic] = numberRowBasic;
389 // can change for gub so redo
390 numberColumnBasic = numberBasic - numberRowBasic;
391 matrix->fillBasis(model,
392 pivotTemp + numberRowBasic,
393 numberColumnBasic,
394 indexRowU,
395 startColumnU + numberRowBasic,
396 numberInRow,
397 numberInColumn + numberRowBasic,
398 elementU);
399#if 0
400 {
401 printf("%d row basic, %d column basic\n", numberRowBasic, numberColumnBasic);
402 for (int i = 0; i < numberElements; i++)
403 printf("row %d col %d value %g\n", indexRowU_.array()[i], indexColumnU_[i],
404 elementU_.array()[i]);
405 }
406#endif
407 // recompute number basic
408 numberBasic = numberRowBasic + numberColumnBasic;
409 if (numberBasic)
410 numberElements = startColumnU[numberBasic-1]
411 + numberInColumn[numberBasic-1];
412 else
413 numberElements = 0;
414 lengthU_ = numberElements;
415 //saveFactorization("dump.d");
416 if (biasLU_ >= 3 || numberRows_ != numberColumns_)
417 preProcess ( 2 );
418 else
419 preProcess ( 3 ); // no row copy
420 factor ( );
421 if (status_ == -99) {
422 // get more memory
423 areaFactor(2.0 * areaFactor());
424 } else if (status_ == -1 && model->numberIterations() == 0 &&
425 denseThreshold_) {
426 // Round again without dense
427 denseThreshold_ = 0;
428 status_ = -99;
429 }
430 }
431 // If we get here status is 0 or -1
432 if (status_ == 0) {
433 // We may need to tamper with order and redo - e.g. network with side
434 int useNumberRows = numberRows;
435 // **** we will also need to add test in dual steepest to do
436 // as we do for network
437 matrix->generalExpanded(model, 12, useNumberRows);
438 const int * permuteBack = permuteBack_.array();
439 const int * back = pivotColumnBack_.array();
440 //int * pivotTemp = pivotColumn_.array();
441 //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp );
442 // Redo pivot order
443 for (i = 0; i < numberRowBasic; i++) {
444 int k = pivotTemp[i];
445 // so rowIsBasic[k] would be permuteBack[back[i]]
446 pivotVariable[permuteBack[back[i]]] = k + numberColumns;
447 }
448 for (; i < useNumberRows; i++) {
449 int k = pivotTemp[i];
450 // so rowIsBasic[k] would be permuteBack[back[i]]
451 pivotVariable[permuteBack[back[i]]] = k;
452 }
453#if 0
454 if (numberSave >= 0) {
455 numberSave = numberDense_;
456 memset(saveList, 0, ((numberRows_ + 31) >> 5)*sizeof(int));
457 for (i = numberRows_ - numberSave; i < numberRows_; i++) {
458 int k = pivotTemp[pivotColumn_.array()[i]];
459 setDense(k);
460 }
461 }
462#endif
463 // Set up permutation vector
464 // these arrays start off as copies of permute
465 // (and we could use permute_ instead of pivotColumn (not back though))
466 ClpDisjointCopyN ( permute_.array(), useNumberRows , pivotColumn_.array() );
467 ClpDisjointCopyN ( permuteBack_.array(), useNumberRows , pivotColumnBack_.array() );
468#ifndef SLIM_CLP
469 if (networkMatrix) {
470 maximumPivots(CoinMax(2000, maximumPivots()));
471 // redo arrays
472 for (int iRow = 0; iRow < 4; iRow++) {
473 int length = model->numberRows() + maximumPivots();
474 if (iRow == 3 || model->objectiveAsObject()->type() > 1)
475 length += model->numberColumns();
476 model->rowArray(iRow)->reserve(length);
477 }
478 // create network factorization
479 if (doCheck)
480 delete networkBasis_; // temp
481 networkBasis_ = new ClpNetworkBasis(model, numberRows_,
482 pivotRegion_.array(),
483 permuteBack_.array(),
484 startColumnU_.array(),
485 numberInColumn_.array(),
486 indexRowU_.array(),
487 elementU_.array());
488 // kill off arrays in ordinary factorization
489 if (!doCheck) {
490 gutsOfDestructor();
491 // but make sure numberRows_ set
492 numberRows_ = model->numberRows();
493 status_ = 0;
494#if 0
495 // but put back permute arrays so odd things will work
496 int numberRows = model->numberRows();
497 pivotColumnBack_ = new int [numberRows];
498 permute_ = new int [numberRows];
499 int i;
500 for (i = 0; i < numberRows; i++) {
501 pivotColumnBack_[i] = i;
502 permute_[i] = i;
503 }
504#endif
505 }
506 } else {
507#endif
508 // See if worth going sparse and when
509 if (numberFtranCounts_ > 100) {
510 ftranCountInput_ = CoinMax(ftranCountInput_, 1.0);
511 ftranAverageAfterL_ = CoinMax(ftranCountAfterL_ / ftranCountInput_, 1.0);
512 ftranAverageAfterR_ = CoinMax(ftranCountAfterR_ / ftranCountAfterL_, 1.0);
513 ftranAverageAfterU_ = CoinMax(ftranCountAfterU_ / ftranCountAfterR_, 1.0);
514 if (btranCountInput_ && btranCountAfterU_ && btranCountAfterR_) {
515 btranAverageAfterU_ = CoinMax(btranCountAfterU_ / btranCountInput_, 1.0);
516 btranAverageAfterR_ = CoinMax(btranCountAfterR_ / btranCountAfterU_, 1.0);
517 btranAverageAfterL_ = CoinMax(btranCountAfterL_ / btranCountAfterR_, 1.0);
518 } else {
519 // we have not done any useful btrans (values pass?)
520 btranAverageAfterU_ = 1.0;
521 btranAverageAfterR_ = 1.0;
522 btranAverageAfterL_ = 1.0;
523 }
524 }
525 // scale back
526
527 ftranCountInput_ *= 0.8;
528 ftranCountAfterL_ *= 0.8;
529 ftranCountAfterR_ *= 0.8;
530 ftranCountAfterU_ *= 0.8;
531 btranCountInput_ *= 0.8;
532 btranCountAfterU_ *= 0.8;
533 btranCountAfterR_ *= 0.8;
534 btranCountAfterL_ *= 0.8;
535#ifndef SLIM_CLP
536 }
537#endif
538 } else if (status_ == -1 && (solveType == 0 || solveType == 2)) {
539 // This needs redoing as it was merged coding - does not need array
540 int numberTotal = numberRows + numberColumns;
541 int * isBasic = new int [numberTotal];
542 int * rowIsBasic = isBasic + numberColumns;
543 int * columnIsBasic = isBasic;
544 for (i = 0; i < numberTotal; i++)
545 isBasic[i] = -1;
546 for (i = 0; i < numberRowBasic; i++) {
547 int iRow = pivotTemp[i];
548 rowIsBasic[iRow] = 1;
549 }
550 for (; i < numberBasic; i++) {
551 int iColumn = pivotTemp[i];
552 columnIsBasic[iColumn] = 1;
553 }
554 numberBasic = 0;
555 for (i = 0; i < numberRows; i++)
556 pivotVariable[i] = -1;
557 // mark as basic or non basic
558 const int * pivotColumn = pivotColumn_.array();
559 for (i = 0; i < numberRows; i++) {
560 if (rowIsBasic[i] >= 0) {
561 if (pivotColumn[numberBasic] >= 0) {
562 rowIsBasic[i] = pivotColumn[numberBasic];
563 } else {
564 rowIsBasic[i] = -1;
565 model->setRowStatus(i, ClpSimplex::superBasic);
566 }
567 numberBasic++;
568 }
569 }
570 for (i = 0; i < numberColumns; i++) {
571 if (columnIsBasic[i] >= 0) {
572 if (pivotColumn[numberBasic] >= 0)
573 columnIsBasic[i] = pivotColumn[numberBasic];
574 else
575 columnIsBasic[i] = -1;
576 numberBasic++;
577 }
578 }
579 // leave pivotVariable in useful form for cleaning basis
580 int * pivotVariable = model->pivotVariable();
581 for (i = 0; i < numberRows; i++) {
582 pivotVariable[i] = -1;
583 }
584
585 for (i = 0; i < numberRows; i++) {
586 if (model->getRowStatus(i) == ClpSimplex::basic) {
587 int iPivot = rowIsBasic[i];
588 if (iPivot >= 0)
589 pivotVariable[iPivot] = i + numberColumns;
590 }
591 }
592 for (i = 0; i < numberColumns; i++) {
593 if (model->getColumnStatus(i) == ClpSimplex::basic) {
594 int iPivot = columnIsBasic[i];
595 if (iPivot >= 0)
596 pivotVariable[iPivot] = i;
597 }
598 }
599 delete [] isBasic;
600 double * columnLower = model->lowerRegion();
601 double * columnUpper = model->upperRegion();
602 double * columnActivity = model->solutionRegion();
603 double * rowLower = model->lowerRegion(0);
604 double * rowUpper = model->upperRegion(0);
605 double * rowActivity = model->solutionRegion(0);
606 //redo basis - first take ALL columns out
607 int iColumn;
608 double largeValue = model->largeValue();
609 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
610 if (model->getColumnStatus(iColumn) == ClpSimplex::basic) {
611 // take out
612 if (!valuesPass) {
613 double lower = columnLower[iColumn];
614 double upper = columnUpper[iColumn];
615 double value = columnActivity[iColumn];
616 if (lower > -largeValue || upper < largeValue) {
617 if (fabs(value - lower) < fabs(value - upper)) {
618 model->setColumnStatus(iColumn, ClpSimplex::atLowerBound);
619 columnActivity[iColumn] = lower;
620 } else {
621 model->setColumnStatus(iColumn, ClpSimplex::atUpperBound);
622 columnActivity[iColumn] = upper;
623 }
624 } else {
625 model->setColumnStatus(iColumn, ClpSimplex::isFree);
626 }
627 } else {
628 model->setColumnStatus(iColumn, ClpSimplex::superBasic);
629 }
630 }
631 }
632 int iRow;
633 for (iRow = 0; iRow < numberRows; iRow++) {
634 int iSequence = pivotVariable[iRow];
635 if (iSequence >= 0) {
636 // basic
637 if (iSequence >= numberColumns) {
638 // slack in - leave
639 //assert (iSequence-numberColumns==iRow);
640 } else {
641 assert(model->getRowStatus(iRow) != ClpSimplex::basic);
642 // put back structural
643 model->setColumnStatus(iSequence, ClpSimplex::basic);
644 }
645 } else {
646 // put in slack
647 model->setRowStatus(iRow, ClpSimplex::basic);
648 }
649 }
650 // Put back any key variables for gub (status_ not touched)
651 matrix->generalExpanded(model, 1, status_);
652 // signal repeat
653 status_ = -99;
654 // set fixed if they are
655 for (iRow = 0; iRow < numberRows; iRow++) {
656 if (model->getRowStatus(iRow) != ClpSimplex::basic ) {
657 if (rowLower[iRow] == rowUpper[iRow]) {
658 rowActivity[iRow] = rowLower[iRow];
659 model->setRowStatus(iRow, ClpSimplex::isFixed);
660 }
661 }
662 }
663 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
664 if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) {
665 if (columnLower[iColumn] == columnUpper[iColumn]) {
666 columnActivity[iColumn] = columnLower[iColumn];
667 model->setColumnStatus(iColumn, ClpSimplex::isFixed);
668 }
669 }
670 }
671 }
672 }
673#ifndef SLIM_CLP
674 } else {
675 // network - fake factorization - do nothing
676 status_ = 0;
677 numberPivots_ = 0;
678 }
679#endif
680#ifndef SLIM_CLP
681 if (!status_) {
682 // take out part if quadratic
683 if (model->algorithm() == 2) {
684 ClpObjective * obj = model->objectiveAsObject();
685#ifndef NDEBUG
686 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
687 assert (quadraticObj);
688#else
689 ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
690#endif
691 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
692 int numberXColumns = quadratic->getNumCols();
693 assert (numberXColumns < numberColumns);
694 int base = numberColumns - numberXColumns;
695 int * which = new int [numberXColumns];
696 int * pivotVariable = model->pivotVariable();
697 int * permute = pivotColumn();
698 int i;
699 int n = 0;
700 for (i = 0; i < numberRows; i++) {
701 int iSj = pivotVariable[i] - base;
702 if (iSj >= 0 && iSj < numberXColumns)
703 which[n++] = permute[i];
704 }
705 if (n)
706 emptyRows(n, which);
707 delete [] which;
708 }
709 }
710#endif
711#ifdef CLP_FACTORIZATION_INSTRUMENT
712 factorization_instrument(2);
713#endif
714 return status_;
715}
716/* Replaces one Column to basis,
717 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
718 If checkBeforeModifying is true will do all accuracy checks
719 before modifying factorization. Whether to set this depends on
720 speed considerations. You could just do this on first iteration
721 after factorization and thereafter re-factorize
722 partial update already in U */
723int
724ClpFactorization::replaceColumn ( const ClpSimplex * model,
725 CoinIndexedVector * regionSparse,
726 CoinIndexedVector * tableauColumn,
727 int pivotRow,
728 double pivotCheck ,
729 bool checkBeforeModifying,
730 double acceptablePivot)
731{
732 int returnCode;
733#ifndef SLIM_CLP
734 if (!networkBasis_) {
735#endif
736#ifdef CLP_FACTORIZATION_INSTRUMENT
737 factorization_instrument(-1);
738#endif
739 // see if FT
740 if (doForrestTomlin_) {
741 returnCode = CoinFactorization::replaceColumn(regionSparse,
742 pivotRow,
743 pivotCheck,
744 checkBeforeModifying,
745 acceptablePivot);
746 } else {
747 returnCode = CoinFactorization::replaceColumnPFI(tableauColumn,
748 pivotRow, pivotCheck); // Note array
749 }
750#ifdef CLP_FACTORIZATION_INSTRUMENT
751 factorization_instrument(3);
752#endif
753
754#ifndef SLIM_CLP
755 } else {
756 if (doCheck) {
757 returnCode = CoinFactorization::replaceColumn(regionSparse,
758 pivotRow,
759 pivotCheck,
760 checkBeforeModifying,
761 acceptablePivot);
762 networkBasis_->replaceColumn(regionSparse,
763 pivotRow);
764 } else {
765 // increase number of pivots
766 numberPivots_++;
767 returnCode = networkBasis_->replaceColumn(regionSparse,
768 pivotRow);
769 }
770 }
771#endif
772 return returnCode;
773}
774
775/* Updates one column (FTRAN) from region2
776 number returned is negative if no room
777 region1 starts as zero and is zero at end */
778int
779ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
780 CoinIndexedVector * regionSparse2)
781{
782#ifdef CLP_DEBUG
783 regionSparse->checkClear();
784#endif
785 if (!numberRows_)
786 return 0;
787#ifndef SLIM_CLP
788 if (!networkBasis_) {
789#endif
790#ifdef CLP_FACTORIZATION_INSTRUMENT
791 factorization_instrument(-1);
792#endif
793 collectStatistics_ = true;
794 int returnCode = CoinFactorization::updateColumnFT(regionSparse,
795 regionSparse2);
796 collectStatistics_ = false;
797#ifdef CLP_FACTORIZATION_INSTRUMENT
798 factorization_instrument(4);
799#endif
800 return returnCode;
801#ifndef SLIM_CLP
802 } else {
803#ifdef CHECK_NETWORK
804 CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
805 double * check = new double[numberRows_];
806 int returnCode = CoinFactorization::updateColumnFT(regionSparse,
807 regionSparse2);
808 networkBasis_->updateColumn(regionSparse, save, -1);
809 int i;
810 double * array = regionSparse2->denseVector();
811 int * indices = regionSparse2->getIndices();
812 int n = regionSparse2->getNumElements();
813 memset(check, 0, numberRows_ * sizeof(double));
814 double * array2 = save->denseVector();
815 int * indices2 = save->getIndices();
816 int n2 = save->getNumElements();
817 assert (n == n2);
818 if (save->packedMode()) {
819 for (i = 0; i < n; i++) {
820 check[indices[i]] = array[i];
821 }
822 for (i = 0; i < n; i++) {
823 double value2 = array2[i];
824 assert (check[indices2[i]] == value2);
825 }
826 } else {
827 for (i = 0; i < numberRows_; i++) {
828 double value1 = array[i];
829 double value2 = array2[i];
830 assert (value1 == value2);
831 }
832 }
833 delete save;
834 delete [] check;
835 return returnCode;
836#else
837 networkBasis_->updateColumn(regionSparse, regionSparse2, -1);
838 return 1;
839#endif
840 }
841#endif
842}
843/* Updates one column (FTRAN) from region2
844 number returned is negative if no room
845 region1 starts as zero and is zero at end */
846int
847ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
848 CoinIndexedVector * regionSparse2,
849 bool noPermute) const
850{
851#ifdef CLP_DEBUG
852 if (!noPermute)
853 regionSparse->checkClear();
854#endif
855 if (!numberRows_)
856 return 0;
857#ifndef SLIM_CLP
858 if (!networkBasis_) {
859#endif
860#ifdef CLP_FACTORIZATION_INSTRUMENT
861 factorization_instrument(-1);
862#endif
863 collectStatistics_ = true;
864 int returnCode = CoinFactorization::updateColumn(regionSparse,
865 regionSparse2,
866 noPermute);
867 collectStatistics_ = false;
868#ifdef CLP_FACTORIZATION_INSTRUMENT
869 factorization_instrument(5);
870#endif
871 return returnCode;
872#ifndef SLIM_CLP
873 } else {
874#ifdef CHECK_NETWORK
875 CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
876 double * check = new double[numberRows_];
877 int returnCode = CoinFactorization::updateColumn(regionSparse,
878 regionSparse2,
879 noPermute);
880 networkBasis_->updateColumn(regionSparse, save, -1);
881 int i;
882 double * array = regionSparse2->denseVector();
883 int * indices = regionSparse2->getIndices();
884 int n = regionSparse2->getNumElements();
885 memset(check, 0, numberRows_ * sizeof(double));
886 double * array2 = save->denseVector();
887 int * indices2 = save->getIndices();
888 int n2 = save->getNumElements();
889 assert (n == n2);
890 if (save->packedMode()) {
891 for (i = 0; i < n; i++) {
892 check[indices[i]] = array[i];
893 }
894 for (i = 0; i < n; i++) {
895 double value2 = array2[i];
896 assert (check[indices2[i]] == value2);
897 }
898 } else {
899 for (i = 0; i < numberRows_; i++) {
900 double value1 = array[i];
901 double value2 = array2[i];
902 assert (value1 == value2);
903 }
904 }
905 delete save;
906 delete [] check;
907 return returnCode;
908#else
909 networkBasis_->updateColumn(regionSparse, regionSparse2, -1);
910 return 1;
911#endif
912 }
913#endif
914}
915/* Updates one column (FTRAN) from region2
916 Tries to do FT update
917 number returned is negative if no room.
918 Also updates region3
919 region1 starts as zero and is zero at end */
920int
921ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
922 CoinIndexedVector * regionSparse2,
923 CoinIndexedVector * regionSparse3,
924 bool noPermuteRegion3)
925{
926 int returnCode = updateColumnFT(regionSparse1, regionSparse2);
927 updateColumn(regionSparse1, regionSparse3, noPermuteRegion3);
928 return returnCode;
929}
930/* Updates one column (FTRAN) from region2
931 number returned is negative if no room
932 region1 starts as zero and is zero at end */
933int
934ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
935 CoinIndexedVector * regionSparse2,
936 bool noPermute) const
937{
938 if (!noPermute)
939 regionSparse->checkClear();
940 if (!numberRows_)
941 return 0;
942 collectStatistics_ = false;
943 int returnCode = CoinFactorization::updateColumn(regionSparse,
944 regionSparse2,
945 noPermute);
946 return returnCode;
947}
948/* Updates one column (BTRAN) from region2
949 region1 starts as zero and is zero at end */
950int
951ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
952 CoinIndexedVector * regionSparse2) const
953{
954 if (!numberRows_)
955 return 0;
956#ifndef SLIM_CLP
957 if (!networkBasis_) {
958#endif
959#ifdef CLP_FACTORIZATION_INSTRUMENT
960 factorization_instrument(-1);
961#endif
962 collectStatistics_ = true;
963 int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
964 regionSparse2);
965 collectStatistics_ = false;
966#ifdef CLP_FACTORIZATION_INSTRUMENT
967 factorization_instrument(6);
968#endif
969 return returnCode;
970#ifndef SLIM_CLP
971 } else {
972#ifdef CHECK_NETWORK
973 CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
974 double * check = new double[numberRows_];
975 int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
976 regionSparse2);
977 networkBasis_->updateColumnTranspose(regionSparse, save);
978 int i;
979 double * array = regionSparse2->denseVector();
980 int * indices = regionSparse2->getIndices();
981 int n = regionSparse2->getNumElements();
982 memset(check, 0, numberRows_ * sizeof(double));
983 double * array2 = save->denseVector();
984 int * indices2 = save->getIndices();
985 int n2 = save->getNumElements();
986 assert (n == n2);
987 if (save->packedMode()) {
988 for (i = 0; i < n; i++) {
989 check[indices[i]] = array[i];
990 }
991 for (i = 0; i < n; i++) {
992 double value2 = array2[i];
993 assert (check[indices2[i]] == value2);
994 }
995 } else {
996 for (i = 0; i < numberRows_; i++) {
997 double value1 = array[i];
998 double value2 = array2[i];
999 assert (value1 == value2);
1000 }
1001 }
1002 delete save;
1003 delete [] check;
1004 return returnCode;
1005#else
1006 return networkBasis_->updateColumnTranspose(regionSparse, regionSparse2);
1007#endif
1008 }
1009#endif
1010}
1011/* makes a row copy of L for speed and to allow very sparse problems */
1012void
1013ClpFactorization::goSparse()
1014{
1015#ifdef CLP_FACTORIZATION_INSTRUMENT
1016 factorization_instrument(-1);
1017#endif
1018#ifndef SLIM_CLP
1019 if (!networkBasis_)
1020#endif
1021 CoinFactorization::goSparse();
1022#ifdef CLP_FACTORIZATION_INSTRUMENT
1023 factorization_instrument(7);
1024#endif
1025}
1026// Cleans up i.e. gets rid of network basis
1027void
1028ClpFactorization::cleanUp()
1029{
1030#ifndef SLIM_CLP
1031 delete networkBasis_;
1032 networkBasis_ = NULL;
1033#endif
1034 resetStatistics();
1035}
1036/// Says whether to redo pivot order
1037bool
1038ClpFactorization::needToReorder() const
1039{
1040#ifdef CHECK_NETWORK
1041 return true;
1042#endif
1043#ifndef SLIM_CLP
1044 if (!networkBasis_)
1045#endif
1046 return true;
1047#ifndef SLIM_CLP
1048 else
1049 return false;
1050#endif
1051}
1052// Get weighted row list
1053void
1054ClpFactorization::getWeights(int * weights) const
1055{
1056#ifdef CLP_FACTORIZATION_INSTRUMENT
1057 factorization_instrument(-1);
1058#endif
1059#ifndef SLIM_CLP
1060 if (networkBasis_) {
1061 // Network - just unit
1062 for (int i = 0; i < numberRows_; i++)
1063 weights[i] = 1;
1064 return;
1065 }
1066#endif
1067 int * numberInRow = numberInRow_.array();
1068 int * numberInColumn = numberInColumn_.array();
1069 int * permuteBack = pivotColumnBack_.array();
1070 int * indexRowU = indexRowU_.array();
1071 const CoinBigIndex * startColumnU = startColumnU_.array();
1072 const CoinBigIndex * startRowL = startRowL_.array();
1073 if (!startRowL || !numberInRow_.array()) {
1074 int * temp = new int[numberRows_];
1075 memset(temp, 0, numberRows_ * sizeof(int));
1076 int i;
1077 for (i = 0; i < numberRows_; i++) {
1078 // one for pivot
1079 temp[i]++;
1080 CoinBigIndex j;
1081 for (j = startColumnU[i]; j < startColumnU[i] + numberInColumn[i]; j++) {
1082 int iRow = indexRowU[j];
1083 temp[iRow]++;
1084 }
1085 }
1086 CoinBigIndex * startColumnL = startColumnL_.array();
1087 int * indexRowL = indexRowL_.array();
1088 for (i = baseL_; i < baseL_ + numberL_; i++) {
1089 CoinBigIndex j;
1090 for (j = startColumnL[i]; j < startColumnL[i+1]; j++) {
1091 int iRow = indexRowL[j];
1092 temp[iRow]++;
1093 }
1094 }
1095 for (i = 0; i < numberRows_; i++) {
1096 int number = temp[i];
1097 int iPermute = permuteBack[i];
1098 weights[iPermute] = number;
1099 }
1100 delete [] temp;
1101 } else {
1102 int i;
1103 for (i = 0; i < numberRows_; i++) {
1104 int number = startRowL[i+1] - startRowL[i] + numberInRow[i] + 1;
1105 //number = startRowL[i+1]-startRowL[i]+1;
1106 //number = numberInRow[i]+1;
1107 int iPermute = permuteBack[i];
1108 weights[iPermute] = number;
1109 }
1110 }
1111#ifdef CLP_FACTORIZATION_INSTRUMENT
1112 factorization_instrument(8);
1113#endif
1114}
1115#else
1116// This one allows multiple factorizations
1117#if CLP_MULTIPLE_FACTORIZATIONS == 1
1118typedef CoinDenseFactorization CoinOtherFactorization;
1119typedef CoinOslFactorization CoinOtherFactorization;
1120#elif CLP_MULTIPLE_FACTORIZATIONS == 2
1121#include "CoinSimpFactorization.hpp"
1122typedef CoinSimpFactorization CoinOtherFactorization;
1123typedef CoinOslFactorization CoinOtherFactorization;
1124#elif CLP_MULTIPLE_FACTORIZATIONS == 3
1125#include "CoinSimpFactorization.hpp"
1126#define CoinOslFactorization CoinDenseFactorization
1127#elif CLP_MULTIPLE_FACTORIZATIONS == 4
1128#include "CoinSimpFactorization.hpp"
1129//#define CoinOslFactorization CoinDenseFactorization
1130#include "CoinOslFactorization.hpp"
1131#endif
1132
1133//-------------------------------------------------------------------
1134// Default Constructor
1135//-------------------------------------------------------------------
1136ClpFactorization::ClpFactorization ()
1137{
1138#ifndef SLIM_CLP
1139 networkBasis_ = NULL;
1140#endif
1141 //coinFactorizationA_ = NULL;
1142 coinFactorizationA_ = new CoinFactorization() ;
1143 coinFactorizationB_ = NULL;
1144 //coinFactorizationB_ = new CoinOtherFactorization();
1145 forceB_ = 0;
1146 goOslThreshold_ = -1;
1147 goDenseThreshold_ = -1;
1148 goSmallThreshold_ = -1;
1149}
1150
1151//-------------------------------------------------------------------
1152// Copy constructor
1153//-------------------------------------------------------------------
1154ClpFactorization::ClpFactorization (const ClpFactorization & rhs,
1155 int denseIfSmaller)
1156{
1157#ifdef CLP_FACTORIZATION_INSTRUMENT
1158 factorization_instrument(-1);
1159#endif
1160#ifndef SLIM_CLP
1161 if (rhs.networkBasis_)
1162 networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1163 else
1164 networkBasis_ = NULL;
1165#endif
1166 forceB_ = rhs.forceB_;
1167 goOslThreshold_ = rhs.goOslThreshold_;
1168 goDenseThreshold_ = rhs.goDenseThreshold_;
1169 goSmallThreshold_ = rhs.goSmallThreshold_;
1170 int goDense = 0;
1171#ifdef CLP_REUSE_ETAS
1172 model_=rhs.model_;
1173#endif
1174 if (denseIfSmaller > 0 && denseIfSmaller <= goDenseThreshold_) {
1175 CoinDenseFactorization * denseR =
1176 dynamic_cast<CoinDenseFactorization *>(rhs.coinFactorizationB_);
1177 if (!denseR)
1178 goDense = 1;
1179 }
1180 if (denseIfSmaller > 0 && !rhs.coinFactorizationB_) {
1181 if (denseIfSmaller <= goDenseThreshold_)
1182 goDense = 1;
1183 else if (denseIfSmaller <= goSmallThreshold_)
1184 goDense = 2;
1185 else if (denseIfSmaller <= goOslThreshold_)
1186 goDense = 3;
1187 } else if (denseIfSmaller < 0) {
1188 if (-denseIfSmaller <= goDenseThreshold_)
1189 goDense = 1;
1190 else if (-denseIfSmaller <= goSmallThreshold_)
1191 goDense = 2;
1192 else if (-denseIfSmaller <= goOslThreshold_)
1193 goDense = 3;
1194 }
1195 if (rhs.coinFactorizationA_ && !goDense)
1196 coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
1197 else
1198 coinFactorizationA_ = NULL;
1199 if (rhs.coinFactorizationB_ && (denseIfSmaller >= 0 || !goDense))
1200 coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1201 else
1202 coinFactorizationB_ = NULL;
1203 if (goDense) {
1204 delete coinFactorizationB_;
1205 if (goDense == 1)
1206 coinFactorizationB_ = new CoinDenseFactorization();
1207 else if (goDense == 2)
1208 coinFactorizationB_ = new CoinSimpFactorization();
1209 else
1210 coinFactorizationB_ = new CoinOslFactorization();
1211 if (rhs.coinFactorizationA_) {
1212 coinFactorizationB_->maximumPivots(rhs.coinFactorizationA_->maximumPivots());
1213 coinFactorizationB_->pivotTolerance(rhs.coinFactorizationA_->pivotTolerance());
1214 coinFactorizationB_->zeroTolerance(rhs.coinFactorizationA_->zeroTolerance());
1215 } else {
1216 assert (coinFactorizationB_);
1217 coinFactorizationB_->maximumPivots(rhs.coinFactorizationB_->maximumPivots());
1218 coinFactorizationB_->pivotTolerance(rhs.coinFactorizationB_->pivotTolerance());
1219 coinFactorizationB_->zeroTolerance(rhs.coinFactorizationB_->zeroTolerance());
1220 }
1221 }
1222 assert (!coinFactorizationA_ || !coinFactorizationB_);
1223#ifdef CLP_FACTORIZATION_INSTRUMENT
1224 factorization_instrument(1);
1225#endif
1226}
1227
1228ClpFactorization::ClpFactorization (const CoinFactorization & rhs)
1229{
1230#ifdef CLP_FACTORIZATION_INSTRUMENT
1231 factorization_instrument(-1);
1232#endif
1233#ifndef SLIM_CLP
1234 networkBasis_ = NULL;
1235#endif
1236 coinFactorizationA_ = new CoinFactorization(rhs);
1237 coinFactorizationB_ = NULL;
1238#ifdef CLP_FACTORIZATION_INSTRUMENT
1239 factorization_instrument(1);
1240#endif
1241 forceB_ = 0;
1242 goOslThreshold_ = -1;
1243 goDenseThreshold_ = -1;
1244 goSmallThreshold_ = -1;
1245 assert (!coinFactorizationA_ || !coinFactorizationB_);
1246}
1247
1248ClpFactorization::ClpFactorization (const CoinOtherFactorization & rhs)
1249{
1250#ifdef CLP_FACTORIZATION_INSTRUMENT
1251 factorization_instrument(-1);
1252#endif
1253#ifndef SLIM_CLP
1254 networkBasis_ = NULL;
1255#endif
1256 coinFactorizationA_ = NULL;
1257 coinFactorizationB_ = rhs.clone();
1258 //coinFactorizationB_ = new CoinOtherFactorization(rhs);
1259 forceB_ = 0;
1260 goOslThreshold_ = -1;
1261 goDenseThreshold_ = -1;
1262 goSmallThreshold_ = -1;
1263#ifdef CLP_FACTORIZATION_INSTRUMENT
1264 factorization_instrument(1);
1265#endif
1266 assert (!coinFactorizationA_ || !coinFactorizationB_);
1267}
1268
1269//-------------------------------------------------------------------
1270// Destructor
1271//-------------------------------------------------------------------
1272ClpFactorization::~ClpFactorization ()
1273{
1274#ifndef SLIM_CLP
1275 delete networkBasis_;
1276#endif
1277 delete coinFactorizationA_;
1278 delete coinFactorizationB_;
1279}
1280
1281//----------------------------------------------------------------
1282// Assignment operator
1283//-------------------------------------------------------------------
1284ClpFactorization &
1285ClpFactorization::operator=(const ClpFactorization& rhs)
1286{
1287#ifdef CLP_FACTORIZATION_INSTRUMENT
1288 factorization_instrument(-1);
1289#endif
1290 if (this != &rhs) {
1291#ifndef SLIM_CLP
1292 delete networkBasis_;
1293 if (rhs.networkBasis_)
1294 networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1295 else
1296 networkBasis_ = NULL;
1297#endif
1298 forceB_ = rhs.forceB_;
1299#ifdef CLP_REUSE_ETAS
1300 model_=rhs.model_;
1301#endif
1302 goOslThreshold_ = rhs.goOslThreshold_;
1303 goDenseThreshold_ = rhs.goDenseThreshold_;
1304 goSmallThreshold_ = rhs.goSmallThreshold_;
1305 if (rhs.coinFactorizationA_) {
1306 if (coinFactorizationA_)
1307 *coinFactorizationA_ = *(rhs.coinFactorizationA_);
1308 else
1309 coinFactorizationA_ = new CoinFactorization(*rhs.coinFactorizationA_);
1310 } else {
1311 delete coinFactorizationA_;
1312 coinFactorizationA_ = NULL;
1313 }
1314
1315 if (rhs.coinFactorizationB_) {
1316 if (coinFactorizationB_) {
1317 CoinDenseFactorization * denseR = dynamic_cast<CoinDenseFactorization *>(rhs.coinFactorizationB_);
1318 CoinDenseFactorization * dense = dynamic_cast<CoinDenseFactorization *>(coinFactorizationB_);
1319 CoinOslFactorization * oslR = dynamic_cast<CoinOslFactorization *>(rhs.coinFactorizationB_);
1320 CoinOslFactorization * osl = dynamic_cast<CoinOslFactorization *>(coinFactorizationB_);
1321 CoinSimpFactorization * simpR = dynamic_cast<CoinSimpFactorization *>(rhs.coinFactorizationB_);
1322 CoinSimpFactorization * simp = dynamic_cast<CoinSimpFactorization *>(coinFactorizationB_);
1323 if (dense && denseR) {
1324 *dense = *denseR;
1325 } else if (osl && oslR) {
1326 *osl = *oslR;
1327 } else if (simp && simpR) {
1328 *simp = *simpR;
1329 } else {
1330 delete coinFactorizationB_;
1331 coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1332 }
1333 } else {
1334 coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1335 }
1336 } else {
1337 delete coinFactorizationB_;
1338 coinFactorizationB_ = NULL;
1339 }
1340 }
1341#ifdef CLP_FACTORIZATION_INSTRUMENT
1342 factorization_instrument(1);
1343#endif
1344 assert (!coinFactorizationA_ || !coinFactorizationB_);
1345 return *this;
1346}
1347// Go over to dense code
1348void
1349ClpFactorization::goDenseOrSmall(int numberRows)
1350{
1351 if (!forceB_) {
1352 if (numberRows <= goDenseThreshold_) {
1353 delete coinFactorizationA_;
1354 delete coinFactorizationB_;
1355 coinFactorizationA_ = NULL;
1356 coinFactorizationB_ = new CoinDenseFactorization();
1357 //printf("going dense\n");
1358 } else if (numberRows <= goSmallThreshold_) {
1359 delete coinFactorizationA_;
1360 delete coinFactorizationB_;
1361 coinFactorizationA_ = NULL;
1362 coinFactorizationB_ = new CoinSimpFactorization();
1363 //printf("going small\n");
1364 } else if (numberRows <= goOslThreshold_) {
1365 delete coinFactorizationA_;
1366 delete coinFactorizationB_;
1367 coinFactorizationA_ = NULL;
1368 coinFactorizationB_ = new CoinOslFactorization();
1369 //printf("going small\n");
1370 }
1371 }
1372 assert (!coinFactorizationA_ || !coinFactorizationB_);
1373}
1374// If nonzero force use of 1,dense 2,small 3,osl
1375void
1376ClpFactorization::forceOtherFactorization(int which)
1377{
1378 delete coinFactorizationB_;
1379 forceB_ = 0;
1380 coinFactorizationB_ = NULL;
1381 if (which > 0 && which < 4) {
1382 delete coinFactorizationA_;
1383 coinFactorizationA_ = NULL;
1384 forceB_ = which;
1385 switch (which) {
1386 case 1:
1387 coinFactorizationB_ = new CoinDenseFactorization();
1388 goDenseThreshold_ = COIN_INT_MAX;
1389 break;
1390 case 2:
1391 coinFactorizationB_ = new CoinSimpFactorization();
1392 goSmallThreshold_ = COIN_INT_MAX;
1393 break;
1394 case 3:
1395 coinFactorizationB_ = new CoinOslFactorization();
1396 goOslThreshold_ = COIN_INT_MAX;
1397 break;
1398 }
1399 } else if (!coinFactorizationA_) {
1400 coinFactorizationA_ = new CoinFactorization();
1401 goOslThreshold_ = -1;
1402 goDenseThreshold_ = -1;
1403 goSmallThreshold_ = -1;
1404 }
1405}
1406int
1407ClpFactorization::factorize ( ClpSimplex * model,
1408 int solveType, bool valuesPass)
1409{
1410#ifdef CLP_REUSE_ETAS
1411 model_= model;
1412#endif
1413 //if ((model->specialOptions()&16384))
1414 //printf("factor at %d iterations\n",model->numberIterations());
1415 ClpMatrixBase * matrix = model->clpMatrix();
1416 int numberRows = model->numberRows();
1417 int numberColumns = model->numberColumns();
1418 if (!numberRows)
1419 return 0;
1420#ifdef CLP_FACTORIZATION_INSTRUMENT
1421 factorization_instrument(-1);
1422#endif
1423 bool anyChanged = false;
1424 if (coinFactorizationB_) {
1425 coinFactorizationB_->setStatus(-99);
1426 int * pivotVariable = model->pivotVariable();
1427 //returns 0 -okay, -1 singular, -2 too many in basis */
1428 // allow dense
1429 int solveMode = 2;
1430 if (model->numberIterations())
1431 solveMode += 8;
1432 if (valuesPass)
1433 solveMode += 4;
1434 coinFactorizationB_->setSolveMode(solveMode);
1435 while (status() < -98) {
1436
1437 int i;
1438 int numberBasic = 0;
1439 int numberRowBasic;
1440 // Move pivot variables across if they look good
1441 int * pivotTemp = model->rowArray(0)->getIndices();
1442 assert (!model->rowArray(0)->getNumElements());
1443 if (!matrix->rhsOffset(model)) {
1444 // Seems to prefer things in order so quickest
1445 // way is to go though like this
1446 for (i = 0; i < numberRows; i++) {
1447 if (model->getRowStatus(i) == ClpSimplex::basic)
1448 pivotTemp[numberBasic++] = i;
1449 }
1450 numberRowBasic = numberBasic;
1451 /* Put column basic variables into pivotVariable
1452 This is done by ClpMatrixBase to allow override for gub
1453 */
1454 matrix->generalExpanded(model, 0, numberBasic);
1455 } else {
1456 // Long matrix - do a different way
1457 bool fullSearch = false;
1458 for (i = 0; i < numberRows; i++) {
1459 int iPivot = pivotVariable[i];
1460 if (iPivot >= numberColumns) {
1461 pivotTemp[numberBasic++] = iPivot - numberColumns;
1462 }
1463 }
1464 numberRowBasic = numberBasic;
1465 for (i = 0; i < numberRows; i++) {
1466 int iPivot = pivotVariable[i];
1467 if (iPivot < numberColumns) {
1468 if (iPivot >= 0) {
1469 pivotTemp[numberBasic++] = iPivot;
1470 } else {
1471 // not full basis
1472 fullSearch = true;
1473 break;
1474 }
1475 }
1476 }
1477 if (fullSearch) {
1478 // do slow way
1479 numberBasic = 0;
1480 for (i = 0; i < numberRows; i++) {
1481 if (model->getRowStatus(i) == ClpSimplex::basic)
1482 pivotTemp[numberBasic++] = i;
1483 }
1484 numberRowBasic = numberBasic;
1485 /* Put column basic variables into pivotVariable
1486 This is done by ClpMatrixBase to allow override for gub
1487 */
1488 matrix->generalExpanded(model, 0, numberBasic);
1489 }
1490 }
1491 if (numberBasic > model->maximumBasic()) {
1492 // Take out some
1493 numberBasic = numberRowBasic;
1494 for (int i = 0; i < numberColumns; i++) {
1495 if (model->getColumnStatus(i) == ClpSimplex::basic) {
1496 if (numberBasic < numberRows)
1497 numberBasic++;
1498 else
1499 model->setColumnStatus(i, ClpSimplex::superBasic);
1500 }
1501 }
1502 numberBasic = numberRowBasic;
1503 matrix->generalExpanded(model, 0, numberBasic);
1504 } else if (numberBasic < numberRows) {
1505 // add in some
1506 int needed = numberRows - numberBasic;
1507 // move up columns
1508 for (i = numberBasic - 1; i >= numberRowBasic; i--)
1509 pivotTemp[i+needed] = pivotTemp[i];
1510 numberRowBasic = 0;
1511 numberBasic = numberRows;
1512 for (i = 0; i < numberRows; i++) {
1513 if (model->getRowStatus(i) == ClpSimplex::basic) {
1514 pivotTemp[numberRowBasic++] = i;
1515 } else if (needed) {
1516 needed--;
1517 model->setRowStatus(i, ClpSimplex::basic);
1518 pivotTemp[numberRowBasic++] = i;
1519 }
1520 }
1521 }
1522 CoinBigIndex numberElements = numberRowBasic;
1523
1524 // compute how much in basis
1525 // can change for gub
1526 int numberColumnBasic = numberBasic - numberRowBasic;
1527
1528 numberElements += matrix->countBasis(pivotTemp + numberRowBasic,
1529 numberColumnBasic);
1530 // Not needed for dense
1531 numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1532 int numberIterations = model->numberIterations();
1533 coinFactorizationB_->setUsefulInformation(&numberIterations, 0);
1534 coinFactorizationB_->getAreas ( numberRows,
1535 numberRowBasic + numberColumnBasic, numberElements,
1536 2 * numberElements );
1537 // Fill in counts so we can skip part of preProcess
1538 // This is NOT needed for dense but would be needed for later versions
1539 CoinFactorizationDouble * elementU;
1540 int * indexRowU;
1541 CoinBigIndex * startColumnU;
1542 int * numberInRow;
1543 int * numberInColumn;
1544 elementU = coinFactorizationB_->elements();
1545 indexRowU = coinFactorizationB_->indices();
1546 startColumnU = coinFactorizationB_->starts();
1547#ifndef COIN_FAST_CODE
1548 double slackValue;
1549 slackValue = coinFactorizationB_->slackValue();
1550#else
1551#define slackValue -1.0
1552#endif
1553 numberInRow = coinFactorizationB_->numberInRow();
1554 numberInColumn = coinFactorizationB_->numberInColumn();
1555 CoinZeroN ( numberInRow, numberRows );
1556 CoinZeroN ( numberInColumn, numberRows );
1557 for (i = 0; i < numberRowBasic; i++) {
1558 int iRow = pivotTemp[i];
1559 // Change pivotTemp to correct sequence
1560 pivotTemp[i] = iRow + numberColumns;
1561 indexRowU[i] = iRow;
1562 startColumnU[i] = i;
1563 elementU[i] = slackValue;
1564 numberInRow[iRow] = 1;
1565 numberInColumn[i] = 1;
1566 }
1567 startColumnU[numberRowBasic] = numberRowBasic;
1568 // can change for gub so redo
1569 numberColumnBasic = numberBasic - numberRowBasic;
1570 matrix->fillBasis(model,
1571 pivotTemp + numberRowBasic,
1572 numberColumnBasic,
1573 indexRowU,
1574 startColumnU + numberRowBasic,
1575 numberInRow,
1576 numberInColumn + numberRowBasic,
1577 elementU);
1578 // recompute number basic
1579 numberBasic = numberRowBasic + numberColumnBasic;
1580 for (i = numberBasic; i < numberRows; i++)
1581 pivotTemp[i] = -1; // mark not there
1582 if (numberBasic)
1583 numberElements = startColumnU[numberBasic-1]
1584 + numberInColumn[numberBasic-1];
1585 else
1586 numberElements = 0;
1587 coinFactorizationB_->preProcess ( );
1588 coinFactorizationB_->factor ( );
1589 if (coinFactorizationB_->status() == -1 &&
1590 (coinFactorizationB_->solveMode() % 3) != 0) {
1591 int solveMode = coinFactorizationB_->solveMode();
1592 solveMode -= solveMode % 3; // so bottom will be 0
1593 coinFactorizationB_->setSolveMode(solveMode);
1594 coinFactorizationB_->setStatus(-99);
1595 }
1596 if (coinFactorizationB_->status() == -99)
1597 continue;
1598 // If we get here status is 0 or -1
1599 if (coinFactorizationB_->status() == 0 && numberBasic == numberRows) {
1600 coinFactorizationB_->postProcess(pivotTemp, pivotVariable);
1601 } else if (solveType == 0 || solveType == 2/*||solveType==1*/) {
1602 // Change pivotTemp to be correct list
1603 anyChanged = true;
1604 coinFactorizationB_->makeNonSingular(pivotTemp, numberColumns);
1605 double * columnLower = model->lowerRegion();
1606 double * columnUpper = model->upperRegion();
1607 double * columnActivity = model->solutionRegion();
1608 double * rowLower = model->lowerRegion(0);
1609 double * rowUpper = model->upperRegion(0);
1610 double * rowActivity = model->solutionRegion(0);
1611 //redo basis - first take ALL out
1612 int iColumn;
1613 double largeValue = model->largeValue();
1614 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1615 if (model->getColumnStatus(iColumn) == ClpSimplex::basic) {
1616 // take out
1617 if (!valuesPass) {
1618 double lower = columnLower[iColumn];
1619 double upper = columnUpper[iColumn];
1620 double value = columnActivity[iColumn];
1621 if (lower > -largeValue || upper < largeValue) {
1622 if (fabs(value - lower) < fabs(value - upper)) {
1623 model->setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1624 columnActivity[iColumn] = lower;
1625 } else {
1626 model->setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1627 columnActivity[iColumn] = upper;
1628 }
1629 } else {
1630 model->setColumnStatus(iColumn, ClpSimplex::isFree);
1631 }
1632 } else {
1633 model->setColumnStatus(iColumn, ClpSimplex::superBasic);
1634 }
1635 }
1636 }
1637 int iRow;
1638 for (iRow = 0; iRow < numberRows; iRow++) {
1639 if (model->getRowStatus(iRow) == ClpSimplex::basic) {
1640 // take out
1641 if (!valuesPass) {
1642 double lower = columnLower[iRow];
1643 double upper = columnUpper[iRow];
1644 double value = columnActivity[iRow];
1645 if (lower > -largeValue || upper < largeValue) {
1646 if (fabs(value - lower) < fabs(value - upper)) {
1647 model->setRowStatus(iRow, ClpSimplex::atLowerBound);
1648 columnActivity[iRow] = lower;
1649 } else {
1650 model->setRowStatus(iRow, ClpSimplex::atUpperBound);
1651 columnActivity[iRow] = upper;
1652 }
1653 } else {
1654 model->setRowStatus(iRow, ClpSimplex::isFree);
1655 }
1656 } else {
1657 model->setRowStatus(iRow, ClpSimplex::superBasic);
1658 }
1659 }
1660 }
1661 for (iRow = 0; iRow < numberRows; iRow++) {
1662 int iSequence = pivotTemp[iRow];
1663 assert (iSequence >= 0);
1664 // basic
1665 model->setColumnStatus(iSequence, ClpSimplex::basic);
1666 }
1667 // signal repeat
1668 coinFactorizationB_->setStatus(-99);
1669 // set fixed if they are
1670 for (iRow = 0; iRow < numberRows; iRow++) {
1671 if (model->getRowStatus(iRow) != ClpSimplex::basic ) {
1672 if (rowLower[iRow] == rowUpper[iRow]) {
1673 rowActivity[iRow] = rowLower[iRow];
1674 model->setRowStatus(iRow, ClpSimplex::isFixed);
1675 }
1676 }
1677 }
1678 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1679 if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) {
1680 if (columnLower[iColumn] == columnUpper[iColumn]) {
1681 columnActivity[iColumn] = columnLower[iColumn];
1682 model->setColumnStatus(iColumn, ClpSimplex::isFixed);
1683 }
1684 }
1685 }
1686 }
1687 }
1688#ifdef CLP_DEBUG
1689 // check basic
1690 CoinIndexedVector region1(2 * numberRows);
1691 CoinIndexedVector region2B(2 * numberRows);
1692 int iPivot;
1693 double * arrayB = region2B.denseVector();
1694 int i;
1695 for (iPivot = 0; iPivot < numberRows; iPivot++) {
1696 int iSequence = pivotVariable[iPivot];
1697 model->unpack(&region2B, iSequence);
1698 coinFactorizationB_->updateColumn(&region1, &region2B);
1699 if (fabs(arrayB[iPivot] - 1.0) < 1.0e-4) {
1700 // OK?
1701 arrayB[iPivot] = 0.0;
1702 } else {
1703 assert (fabs(arrayB[iPivot]) < 1.0e-4);
1704 for (i = 0; i < numberRows; i++) {
1705 if (fabs(arrayB[i] - 1.0) < 1.0e-4)
1706 break;
1707 }
1708 assert (i < numberRows);
1709 printf("variable on row %d landed up on row %d\n", iPivot, i);
1710 arrayB[i] = 0.0;
1711 }
1712 for (i = 0; i < numberRows; i++)
1713 assert (fabs(arrayB[i]) < 1.0e-4);
1714 region2B.clear();
1715 }
1716#endif
1717#ifdef CLP_FACTORIZATION_INSTRUMENT
1718 factorization_instrument(2);
1719#endif
1720 if ( anyChanged && model->algorithm() < 0 && solveType > 0) {
1721 double dummyCost;
1722 static_cast<ClpSimplexDual *> (model)->changeBounds(3,
1723 NULL, dummyCost);
1724 }
1725 return coinFactorizationB_->status();
1726 }
1727 // If too many compressions increase area
1728 if (coinFactorizationA_->pivots() > 1 && coinFactorizationA_->numberCompressions() * 10 > coinFactorizationA_->pivots() + 10) {
1729 coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor() * 1.1);
1730 }
1731 //int numberPivots=coinFactorizationA_->pivots();
1732#if 0
1733 if (model->algorithm() > 0)
1734 numberSave = -1;
1735#endif
1736#ifndef SLIM_CLP
1737 if (!networkBasis_ || doCheck) {
1738#endif
1739 coinFactorizationA_->setStatus(-99);
1740 int * pivotVariable = model->pivotVariable();
1741 int nTimesRound = 0;
1742 //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
1743 while (coinFactorizationA_->status() < -98) {
1744 nTimesRound++;
1745
1746 int i;
1747 int numberBasic = 0;
1748 int numberRowBasic;
1749 // Move pivot variables across if they look good
1750 int * pivotTemp = model->rowArray(0)->getIndices();
1751 assert (!model->rowArray(0)->getNumElements());
1752 if (!matrix->rhsOffset(model)) {
1753#if 0
1754 if (numberSave > 0) {
1755 int nStill = 0;
1756 int nAtBound = 0;
1757 int nZeroDual = 0;
1758 CoinIndexedVector * array = model->rowArray(3);
1759 CoinIndexedVector * objArray = model->columnArray(1);
1760 array->clear();
1761 objArray->clear();
1762 double * cost = model->costRegion();
1763 double tolerance = model->primalTolerance();
1764 double offset = 0.0;
1765 for (i = 0; i < numberRows; i++) {
1766 int iPivot = pivotVariable[i];
1767 if (iPivot < numberColumns && isDense(iPivot)) {
1768 if (model->getColumnStatus(iPivot) == ClpSimplex::basic) {
1769 nStill++;
1770 double value = model->solutionRegion()[iPivot];
1771 double dual = model->dualRowSolution()[i];
1772 double lower = model->lowerRegion()[iPivot];
1773 double upper = model->upperRegion()[iPivot];
1774 ClpSimplex::Status status;
1775 if (fabs(value - lower) < tolerance) {
1776 status = ClpSimplex::atLowerBound;
1777 nAtBound++;
1778 } else if (fabs(value - upper) < tolerance) {
1779 nAtBound++;
1780 status = ClpSimplex::atUpperBound;
1781 } else if (value > lower && value < upper) {
1782 status = ClpSimplex::superBasic;
1783 } else {
1784 status = ClpSimplex::basic;
1785 }
1786 if (status != ClpSimplex::basic) {
1787 if (model->getRowStatus(i) != ClpSimplex::basic) {
1788 model->setColumnStatus(iPivot, ClpSimplex::atLowerBound);
1789 model->setRowStatus(i, ClpSimplex::basic);
1790 pivotVariable[i] = i + numberColumns;
1791 model->dualRowSolution()[i] = 0.0;
1792 model->djRegion(0)[i] = 0.0;
1793 array->add(i, dual);
1794 offset += dual * model->solutionRegion(0)[i];
1795 }
1796 }
1797 if (fabs(dual) < 1.0e-5)
1798 nZeroDual++;
1799 }
1800 }
1801 }
1802 printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
1803 numberSave, nStill, nAtBound, nZeroDual, offset);
1804 if (array->getNumElements()) {
1805 // modify costs
1806 model->clpMatrix()->transposeTimes(model, 1.0, array, model->columnArray(0),
1807 objArray);
1808 array->clear();
1809 int n = objArray->getNumElements();
1810 int * indices = objArray->getIndices();
1811 double * elements = objArray->denseVector();
1812 for (i = 0; i < n; i++) {
1813 int iColumn = indices[i];
1814 cost[iColumn] -= elements[iColumn];
1815 elements[iColumn] = 0.0;
1816 }
1817 objArray->setNumElements(0);
1818 }
1819 }
1820#endif
1821 // Seems to prefer things in order so quickest
1822 // way is to go though like this
1823 for (i = 0; i < numberRows; i++) {
1824 if (model->getRowStatus(i) == ClpSimplex::basic)
1825 pivotTemp[numberBasic++] = i;
1826 }
1827 numberRowBasic = numberBasic;
1828 /* Put column basic variables into pivotVariable
1829 This is done by ClpMatrixBase to allow override for gub
1830 */
1831 matrix->generalExpanded(model, 0, numberBasic);
1832 } else {
1833 // Long matrix - do a different way
1834 bool fullSearch = false;
1835 for (i = 0; i < numberRows; i++) {
1836 int iPivot = pivotVariable[i];
1837 if (iPivot >= numberColumns) {
1838 pivotTemp[numberBasic++] = iPivot - numberColumns;
1839 }
1840 }
1841 numberRowBasic = numberBasic;
1842 for (i = 0; i < numberRows; i++) {
1843 int iPivot = pivotVariable[i];
1844 if (iPivot < numberColumns) {
1845 if (iPivot >= 0) {
1846 pivotTemp[numberBasic++] = iPivot;
1847 } else {
1848 // not full basis
1849 fullSearch = true;
1850 break;
1851 }
1852 }
1853 }
1854 if (fullSearch) {
1855 // do slow way
1856 numberBasic = 0;
1857 for (i = 0; i < numberRows; i++) {
1858 if (model->getRowStatus(i) == ClpSimplex::basic)
1859 pivotTemp[numberBasic++] = i;
1860 }
1861 numberRowBasic = numberBasic;
1862 /* Put column basic variables into pivotVariable
1863 This is done by ClpMatrixBase to allow override for gub
1864 */
1865 matrix->generalExpanded(model, 0, numberBasic);
1866 }
1867 }
1868 if (numberBasic > model->maximumBasic()) {
1869#if 0 // ndef NDEBUG
1870 printf("%d basic - should only be %d\n",
1871 numberBasic, numberRows);
1872#endif
1873 // Take out some
1874 numberBasic = numberRowBasic;
1875 for (int i = 0; i < numberColumns; i++) {
1876 if (model->getColumnStatus(i) == ClpSimplex::basic) {
1877 if (numberBasic < numberRows)
1878 numberBasic++;
1879 else
1880 model->setColumnStatus(i, ClpSimplex::superBasic);
1881 }
1882 }
1883 numberBasic = numberRowBasic;
1884 matrix->generalExpanded(model, 0, numberBasic);
1885 }
1886#ifndef SLIM_CLP
1887 // see if matrix a network
1888#ifndef NO_RTTI
1889 ClpNetworkMatrix* networkMatrix =
1890 dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
1891#else
1892ClpNetworkMatrix* networkMatrix = NULL;
1893if (model->clpMatrix()->type() == 11)
1894 networkMatrix =
1895 static_cast< ClpNetworkMatrix*>(model->clpMatrix());
1896#endif
1897 // If network - still allow ordinary factorization first time for laziness
1898 if (networkMatrix)
1899 coinFactorizationA_->setBiasLU(0); // All to U if network
1900 //int saveMaximumPivots = maximumPivots();
1901 delete networkBasis_;
1902 networkBasis_ = NULL;
1903 if (networkMatrix && !doCheck)
1904 maximumPivots(1);
1905#endif
1906 //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
1907 while (coinFactorizationA_->status() == -99) {
1908 // maybe for speed will be better to leave as many regions as possible
1909 coinFactorizationA_->gutsOfDestructor();
1910 coinFactorizationA_->gutsOfInitialize(2);
1911 CoinBigIndex numberElements = numberRowBasic;
1912
1913 // compute how much in basis
1914
1915 int i;
1916 // can change for gub
1917 int numberColumnBasic = numberBasic - numberRowBasic;
1918
1919 numberElements += matrix->countBasis( pivotTemp + numberRowBasic,
1920 numberColumnBasic);
1921 // and recompute as network side say different
1922 if (model->numberIterations())
1923 numberRowBasic = numberBasic - numberColumnBasic;
1924 numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1925 coinFactorizationA_->getAreas ( numberRows,
1926 numberRowBasic + numberColumnBasic, numberElements,
1927 2 * numberElements );
1928 //fill
1929 // Fill in counts so we can skip part of preProcess
1930 int * numberInRow = coinFactorizationA_->numberInRow();
1931 int * numberInColumn = coinFactorizationA_->numberInColumn();
1932 CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 );
1933 CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 );
1934 CoinFactorizationDouble * elementU = coinFactorizationA_->elementU();
1935 int * indexRowU = coinFactorizationA_->indexRowU();
1936 CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
1937#ifndef COIN_FAST_CODE
1938 double slackValue = coinFactorizationA_->slackValue();
1939#endif
1940 for (i = 0; i < numberRowBasic; i++) {
1941 int iRow = pivotTemp[i];
1942 indexRowU[i] = iRow;
1943 startColumnU[i] = i;
1944 elementU[i] = slackValue;
1945 numberInRow[iRow] = 1;
1946 numberInColumn[i] = 1;
1947 }
1948 startColumnU[numberRowBasic] = numberRowBasic;
1949 // can change for gub so redo
1950 numberColumnBasic = numberBasic - numberRowBasic;
1951 matrix->fillBasis(model,
1952 pivotTemp + numberRowBasic,
1953 numberColumnBasic,
1954 indexRowU,
1955 startColumnU + numberRowBasic,
1956 numberInRow,
1957 numberInColumn + numberRowBasic,
1958 elementU);
1959#if 0
1960 {
1961 printf("%d row basic, %d column basic\n", numberRowBasic, numberColumnBasic);
1962 for (int i = 0; i < numberElements; i++)
1963 printf("row %d col %d value %g\n", indexRowU[i], indexColumnU_[i],
1964 elementU[i]);
1965 }
1966#endif
1967 // recompute number basic
1968 numberBasic = numberRowBasic + numberColumnBasic;
1969 if (numberBasic)
1970 numberElements = startColumnU[numberBasic-1]
1971 + numberInColumn[numberBasic-1];
1972 else
1973 numberElements = 0;
1974 coinFactorizationA_->setNumberElementsU(numberElements);
1975 //saveFactorization("dump.d");
1976 if (coinFactorizationA_->biasLU() >= 3 || coinFactorizationA_->numberRows() != coinFactorizationA_->numberColumns())
1977 coinFactorizationA_->preProcess ( 2 );
1978 else
1979 coinFactorizationA_->preProcess ( 3 ); // no row copy
1980 coinFactorizationA_->factor ( );
1981 if (coinFactorizationA_->status() == -99) {
1982 // get more memory
1983 coinFactorizationA_->areaFactor(2.0 * coinFactorizationA_->areaFactor());
1984 } else if (coinFactorizationA_->status() == -1 &&
1985 (model->numberIterations() == 0 || nTimesRound > 2) &&
1986 coinFactorizationA_->denseThreshold()) {
1987 // Round again without dense
1988 coinFactorizationA_->setDenseThreshold(0);
1989 coinFactorizationA_->setStatus(-99);
1990 }
1991 }
1992 // If we get here status is 0 or -1
1993 if (coinFactorizationA_->status() == 0) {
1994 // We may need to tamper with order and redo - e.g. network with side
1995 int useNumberRows = numberRows;
1996 // **** we will also need to add test in dual steepest to do
1997 // as we do for network
1998 matrix->generalExpanded(model, 12, useNumberRows);
1999 const int * permuteBack = coinFactorizationA_->permuteBack();
2000 const int * back = coinFactorizationA_->pivotColumnBack();
2001 //int * pivotTemp = pivotColumn_.array();
2002 //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp );
2003#ifndef NDEBUG
2004 CoinFillN(pivotVariable, numberRows, -1);
2005#endif
2006 // Redo pivot order
2007 for (i = 0; i < numberRowBasic; i++) {
2008 int k = pivotTemp[i];
2009 // so rowIsBasic[k] would be permuteBack[back[i]]
2010 int j = permuteBack[back[i]];
2011 assert (pivotVariable[j] == -1);
2012 pivotVariable[j] = k + numberColumns;
2013 }
2014 for (; i < useNumberRows; i++) {
2015 int k = pivotTemp[i];
2016 // so rowIsBasic[k] would be permuteBack[back[i]]
2017 int j = permuteBack[back[i]];
2018 assert (pivotVariable[j] == -1);
2019 pivotVariable[j] = k;
2020 }
2021#if 0
2022 if (numberSave >= 0) {
2023 numberSave = numberDense_;
2024 memset(saveList, 0, ((coinFactorizationA_->numberRows() + 31) >> 5)*sizeof(int));
2025 for (i = coinFactorizationA_->numberRows() - numberSave; i < coinFactorizationA_->numberRows(); i++) {
2026 int k = pivotTemp[pivotColumn_.array()[i]];
2027 setDense(k);
2028 }
2029 }
2030#endif
2031 // Set up permutation vector
2032 // these arrays start off as copies of permute
2033 // (and we could use permute_ instead of pivotColumn (not back though))
2034 ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn() );
2035 ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack() );
2036#ifndef SLIM_CLP
2037 if (networkMatrix) {
2038 maximumPivots(CoinMax(2000, maximumPivots()));
2039 // redo arrays
2040 for (int iRow = 0; iRow < 4; iRow++) {
2041 int length = model->numberRows() + maximumPivots();
2042 if (iRow == 3 || model->objectiveAsObject()->type() > 1)
2043 length += model->numberColumns();
2044 model->rowArray(iRow)->reserve(length);
2045 }
2046 // create network factorization
2047 if (doCheck)
2048 delete networkBasis_; // temp
2049 networkBasis_ = new ClpNetworkBasis(model, coinFactorizationA_->numberRows(),
2050 coinFactorizationA_->pivotRegion(),
2051 coinFactorizationA_->permuteBack(),
2052 coinFactorizationA_->startColumnU(),
2053 coinFactorizationA_->numberInColumn(),
2054 coinFactorizationA_->indexRowU(),
2055 coinFactorizationA_->elementU());
2056 // kill off arrays in ordinary factorization
2057 if (!doCheck) {
2058 coinFactorizationA_->gutsOfDestructor();
2059 // but make sure coinFactorizationA_->numberRows() set
2060 coinFactorizationA_->setNumberRows(model->numberRows());
2061 coinFactorizationA_->setStatus(0);
2062#if 0
2063 // but put back permute arrays so odd things will work
2064 int numberRows = model->numberRows();
2065 pivotColumnBack_ = new int [numberRows];
2066 permute_ = new int [numberRows];
2067 int i;
2068 for (i = 0; i < numberRows; i++) {
2069 pivotColumnBack_[i] = i;
2070 permute_[i] = i;
2071 }
2072#endif
2073 }
2074 } else {
2075#endif
2076 // See if worth going sparse and when
2077 coinFactorizationA_->checkSparse();
2078#ifndef SLIM_CLP
2079 }
2080#endif
2081 } else if (coinFactorizationA_->status() == -1 && (solveType == 0 || solveType == 2)) {
2082 // This needs redoing as it was merged coding - does not need array
2083 int numberTotal = numberRows + numberColumns;
2084 int * isBasic = new int [numberTotal];
2085 int * rowIsBasic = isBasic + numberColumns;
2086 int * columnIsBasic = isBasic;
2087 for (i = 0; i < numberTotal; i++)
2088 isBasic[i] = -1;
2089 for (i = 0; i < numberRowBasic; i++) {
2090 int iRow = pivotTemp[i];
2091 rowIsBasic[iRow] = 1;
2092 }
2093 for (; i < numberBasic; i++) {
2094 int iColumn = pivotTemp[i];
2095 columnIsBasic[iColumn] = 1;
2096 }
2097 numberBasic = 0;
2098 for (i = 0; i < numberRows; i++)
2099 pivotVariable[i] = -1;
2100 // mark as basic or non basic
2101 const int * pivotColumn = coinFactorizationA_->pivotColumn();
2102 for (i = 0; i < numberRows; i++) {
2103 if (rowIsBasic[i] >= 0) {
2104 if (pivotColumn[numberBasic] >= 0) {
2105 rowIsBasic[i] = pivotColumn[numberBasic];
2106 } else {
2107 rowIsBasic[i] = -1;
2108 model->setRowStatus(i, ClpSimplex::superBasic);
2109 }
2110 numberBasic++;
2111 }
2112 }
2113 for (i = 0; i < numberColumns; i++) {
2114 if (columnIsBasic[i] >= 0) {
2115 if (pivotColumn[numberBasic] >= 0)
2116 columnIsBasic[i] = pivotColumn[numberBasic];
2117 else
2118 columnIsBasic[i] = -1;
2119 numberBasic++;
2120 }
2121 }
2122 // leave pivotVariable in useful form for cleaning basis
2123 int * pivotVariable = model->pivotVariable();
2124 for (i = 0; i < numberRows; i++) {
2125 pivotVariable[i] = -1;
2126 }
2127
2128 for (i = 0; i < numberRows; i++) {
2129 if (model->getRowStatus(i) == ClpSimplex::basic) {
2130 int iPivot = rowIsBasic[i];
2131 if (iPivot >= 0)
2132 pivotVariable[iPivot] = i + numberColumns;
2133 }
2134 }
2135 for (i = 0; i < numberColumns; i++) {
2136 if (model->getColumnStatus(i) == ClpSimplex::basic) {
2137 int iPivot = columnIsBasic[i];
2138 if (iPivot >= 0)
2139 pivotVariable[iPivot] = i;
2140 }
2141 }
2142 delete [] isBasic;
2143 double * columnLower = model->lowerRegion();
2144 double * columnUpper = model->upperRegion();
2145 double * columnActivity = model->solutionRegion();
2146 double * rowLower = model->lowerRegion(0);
2147 double * rowUpper = model->upperRegion(0);
2148 double * rowActivity = model->solutionRegion(0);
2149 //redo basis - first take ALL columns out
2150 int iColumn;
2151 double largeValue = model->largeValue();
2152 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2153 if (model->getColumnStatus(iColumn) == ClpSimplex::basic) {
2154 // take out
2155 if (!valuesPass) {
2156 double lower = columnLower[iColumn];
2157 double upper = columnUpper[iColumn];
2158 double value = columnActivity[iColumn];
2159 if (lower > -largeValue || upper < largeValue) {
2160 if (fabs(value - lower) < fabs(value - upper)) {
2161 model->setColumnStatus(iColumn, ClpSimplex::atLowerBound);
2162 columnActivity[iColumn] = lower;
2163 } else {
2164 model->setColumnStatus(iColumn, ClpSimplex::atUpperBound);
2165 columnActivity[iColumn] = upper;
2166 }
2167 } else {
2168 model->setColumnStatus(iColumn, ClpSimplex::isFree);
2169 }
2170 } else {
2171 model->setColumnStatus(iColumn, ClpSimplex::superBasic);
2172 }
2173 }
2174 }
2175 int iRow;
2176 for (iRow = 0; iRow < numberRows; iRow++) {
2177 int iSequence = pivotVariable[iRow];
2178 if (iSequence >= 0) {
2179 // basic
2180 if (iSequence >= numberColumns) {
2181 // slack in - leave
2182 //assert (iSequence-numberColumns==iRow);
2183 } else {
2184 assert(model->getRowStatus(iRow) != ClpSimplex::basic);
2185 // put back structural
2186 model->setColumnStatus(iSequence, ClpSimplex::basic);
2187 }
2188 } else {
2189 // put in slack
2190 model->setRowStatus(iRow, ClpSimplex::basic);
2191 }
2192 }
2193 // Put back any key variables for gub
2194 int dummy;
2195 matrix->generalExpanded(model, 1, dummy);
2196 // signal repeat
2197 coinFactorizationA_->setStatus(-99);
2198 // set fixed if they are
2199 for (iRow = 0; iRow < numberRows; iRow++) {
2200 if (model->getRowStatus(iRow) != ClpSimplex::basic ) {
2201 if (rowLower[iRow] == rowUpper[iRow]) {
2202 rowActivity[iRow] = rowLower[iRow];
2203 model->setRowStatus(iRow, ClpSimplex::isFixed);
2204 }
2205 }
2206 }
2207 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2208 if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) {
2209 if (columnLower[iColumn] == columnUpper[iColumn]) {
2210 columnActivity[iColumn] = columnLower[iColumn];
2211 model->setColumnStatus(iColumn, ClpSimplex::isFixed);
2212 }
2213 }
2214 }
2215 }
2216 }
2217#ifndef SLIM_CLP
2218 } else {
2219 // network - fake factorization - do nothing
2220 coinFactorizationA_->setStatus(0);
2221 coinFactorizationA_->setPivots(0);
2222 }
2223#endif
2224#ifndef SLIM_CLP
2225 if (!coinFactorizationA_->status()) {
2226 // take out part if quadratic
2227 if (model->algorithm() == 2) {
2228 ClpObjective * obj = model->objectiveAsObject();
2229#ifndef NDEBUG
2230 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
2231 assert (quadraticObj);
2232#else
2233ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
2234#endif
2235 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2236 int numberXColumns = quadratic->getNumCols();
2237 assert (numberXColumns < numberColumns);
2238 int base = numberColumns - numberXColumns;
2239 int * which = new int [numberXColumns];
2240 int * pivotVariable = model->pivotVariable();
2241 int * permute = pivotColumn();
2242 int i;
2243 int n = 0;
2244 for (i = 0; i < numberRows; i++) {
2245 int iSj = pivotVariable[i] - base;
2246 if (iSj >= 0 && iSj < numberXColumns)
2247 which[n++] = permute[i];
2248 }
2249 if (n)
2250 coinFactorizationA_->emptyRows(n, which);
2251 delete [] which;
2252 }
2253 }
2254#endif
2255#ifdef CLP_FACTORIZATION_INSTRUMENT
2256 factorization_instrument(2);
2257#endif
2258 return coinFactorizationA_->status();
2259}
2260/* Replaces one Column in basis,
2261 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
2262 If checkBeforeModifying is true will do all accuracy checks
2263 before modifying factorization. Whether to set this depends on
2264 speed considerations. You could just do this on first iteration
2265 after factorization and thereafter re-factorize
2266 partial update already in U */
2267int
2268ClpFactorization::replaceColumn ( const ClpSimplex * model,
2269 CoinIndexedVector * regionSparse,
2270 CoinIndexedVector * tableauColumn,
2271 int pivotRow,
2272 double pivotCheck ,
2273 bool checkBeforeModifying,
2274 double acceptablePivot)
2275{
2276#ifndef SLIM_CLP
2277 if (!networkBasis_) {
2278#endif
2279#ifdef CLP_FACTORIZATION_INSTRUMENT
2280 factorization_instrument(-1);
2281#endif
2282 int returnCode;
2283 // see if FT
2284 if (!coinFactorizationA_ || coinFactorizationA_->forrestTomlin()) {
2285 if (coinFactorizationA_) {
2286 returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2287 pivotRow,
2288 pivotCheck,
2289 checkBeforeModifying,
2290 acceptablePivot);
2291 } else {
2292 bool tab = coinFactorizationB_->wantsTableauColumn();
2293#ifdef CLP_REUSE_ETAS
2294 int tempInfo[2];
2295 tempInfo[1] = model_->sequenceOut();
2296#else
2297 int tempInfo[1];
2298#endif
2299 tempInfo[0] = model->numberIterations();
2300 coinFactorizationB_->setUsefulInformation(tempInfo, 1);
2301 returnCode =
2302 coinFactorizationB_->replaceColumn(tab ? tableauColumn : regionSparse,
2303 pivotRow,
2304 pivotCheck,
2305 checkBeforeModifying,
2306 acceptablePivot);
2307#ifdef CLP_DEBUG
2308 // check basic
2309 int numberRows = coinFactorizationB_->numberRows();
2310 CoinIndexedVector region1(2 * numberRows);
2311 CoinIndexedVector region2A(2 * numberRows);
2312 CoinIndexedVector region2B(2 * numberRows);
2313 int iPivot;
2314 double * arrayB = region2B.denseVector();
2315 int * pivotVariable = model->pivotVariable();
2316 int i;
2317 for (iPivot = 0; iPivot < numberRows; iPivot++) {
2318 int iSequence = pivotVariable[iPivot];
2319 if (iPivot == pivotRow)
2320 iSequence = model->sequenceIn();
2321 model->unpack(&region2B, iSequence);
2322 coinFactorizationB_->updateColumn(&region1, &region2B);
2323 assert (fabs(arrayB[iPivot] - 1.0) < 1.0e-4);
2324 arrayB[iPivot] = 0.0;
2325 for (i = 0; i < numberRows; i++)
2326 assert (fabs(arrayB[i]) < 1.0e-4);
2327 region2B.clear();
2328 }
2329#endif
2330 }
2331 } else {
2332 returnCode = coinFactorizationA_->replaceColumnPFI(tableauColumn,
2333 pivotRow, pivotCheck); // Note array
2334 }
2335#ifdef CLP_FACTORIZATION_INSTRUMENT
2336 factorization_instrument(3);
2337#endif
2338 return returnCode;
2339
2340#ifndef SLIM_CLP
2341 } else {
2342 if (doCheck) {
2343 int returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2344 pivotRow,
2345 pivotCheck,
2346 checkBeforeModifying,
2347 acceptablePivot);
2348 networkBasis_->replaceColumn(regionSparse,
2349 pivotRow);
2350 return returnCode;
2351 } else {
2352 // increase number of pivots
2353 coinFactorizationA_->setPivots(coinFactorizationA_->pivots() + 1);
2354 return networkBasis_->replaceColumn(regionSparse,
2355 pivotRow);
2356 }
2357 }
2358#endif
2359}
2360
2361/* Updates one column (FTRAN) from region2
2362 number returned is negative if no room
2363 region1 starts as zero and is zero at end */
2364int
2365ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
2366 CoinIndexedVector * regionSparse2)
2367{
2368#ifdef CLP_DEBUG
2369 regionSparse->checkClear();
2370#endif
2371 if (!numberRows())
2372 return 0;
2373#ifndef SLIM_CLP
2374 if (!networkBasis_) {
2375#endif
2376#ifdef CLP_FACTORIZATION_INSTRUMENT
2377 factorization_instrument(-1);
2378#endif
2379 int returnCode;
2380 if (coinFactorizationA_) {
2381 coinFactorizationA_->setCollectStatistics(true);
2382 returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
2383 regionSparse2);
2384 coinFactorizationA_->setCollectStatistics(false);
2385 } else {
2386#ifdef CLP_REUSE_ETAS
2387 int tempInfo[2];
2388 tempInfo[0] = model_->numberIterations();
2389 tempInfo[1] = model_->sequenceIn();
2390 coinFactorizationB_->setUsefulInformation(tempInfo, 2);
2391#endif
2392 returnCode = coinFactorizationB_->updateColumnFT(regionSparse,
2393 regionSparse2);
2394 }
2395#ifdef CLP_FACTORIZATION_INSTRUMENT
2396 factorization_instrument(4);
2397#endif
2398 return returnCode;
2399#ifndef SLIM_CLP
2400 } else {
2401#ifdef CHECK_NETWORK
2402 CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2403 double * check = new double[coinFactorizationA_->numberRows()];
2404 int returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
2405 regionSparse2);
2406 networkBasis_->updateColumn(regionSparse, save, -1);
2407 int i;
2408 double * array = regionSparse2->denseVector();
2409 int * indices = regionSparse2->getIndices();
2410 int n = regionSparse2->getNumElements();
2411 memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double));
2412 double * array2 = save->denseVector();
2413 int * indices2 = save->getIndices();
2414 int n2 = save->getNumElements();
2415 assert (n == n2);
2416 if (save->packedMode()) {
2417 for (i = 0; i < n; i++) {
2418 check[indices[i]] = array[i];
2419 }
2420 for (i = 0; i < n; i++) {
2421 double value2 = array2[i];
2422 assert (check[indices2[i]] == value2);
2423 }
2424 } else {
2425 int numberRows = coinFactorizationA_->numberRows();
2426 for (i = 0; i < numberRows; i++) {
2427 double value1 = array[i];
2428 double value2 = array2[i];
2429 assert (value1 == value2);
2430 }
2431 }
2432 delete save;
2433 delete [] check;
2434 return returnCode;
2435#else
2436networkBasis_->updateColumn(regionSparse, regionSparse2, -1);
2437return 1;
2438#endif
2439 }
2440#endif
2441}
2442/* Updates one column (FTRAN) from region2
2443 number returned is negative if no room
2444 region1 starts as zero and is zero at end */
2445int
2446ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
2447 CoinIndexedVector * regionSparse2,
2448 bool noPermute) const
2449{
2450#ifdef CLP_DEBUG
2451 if (!noPermute)
2452 regionSparse->checkClear();
2453#endif
2454 if (!numberRows())
2455 return 0;
2456#ifndef SLIM_CLP
2457 if (!networkBasis_) {
2458#endif
2459#ifdef CLP_FACTORIZATION_INSTRUMENT
2460 factorization_instrument(-1);
2461#endif
2462 int returnCode;
2463 if (coinFactorizationA_) {
2464 coinFactorizationA_->setCollectStatistics(true);
2465 returnCode = coinFactorizationA_->updateColumn(regionSparse,
2466 regionSparse2,
2467 noPermute);
2468 coinFactorizationA_->setCollectStatistics(false);
2469 } else {
2470 returnCode = coinFactorizationB_->updateColumn(regionSparse,
2471 regionSparse2,
2472 noPermute);
2473 }
2474#ifdef CLP_FACTORIZATION_INSTRUMENT
2475 factorization_instrument(5);
2476#endif
2477 //#define PRINT_VECTOR
2478#ifdef PRINT_VECTOR
2479 printf("Update\n");
2480 regionSparse2->print();
2481#endif
2482 return returnCode;
2483#ifndef SLIM_CLP
2484 } else {
2485#ifdef CHECK_NETWORK
2486 CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2487 double * check = new double[coinFactorizationA_->numberRows()];
2488 int returnCode = coinFactorizationA_->updateColumn(regionSparse,
2489 regionSparse2,
2490 noPermute);
2491 networkBasis_->updateColumn(regionSparse, save, -1);
2492 int i;
2493 double * array = regionSparse2->denseVector();
2494 int * indices = regionSparse2->getIndices();
2495 int n = regionSparse2->getNumElements();
2496 memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double));
2497 double * array2 = save->denseVector();
2498 int * indices2 = save->getIndices();
2499 int n2 = save->getNumElements();
2500 assert (n == n2);
2501 if (save->packedMode()) {
2502 for (i = 0; i < n; i++) {
2503 check[indices[i]] = array[i];
2504 }
2505 for (i = 0; i < n; i++) {
2506 double value2 = array2[i];
2507 assert (check[indices2[i]] == value2);
2508 }
2509 } else {
2510 int numberRows = coinFactorizationA_->numberRows();
2511 for (i = 0; i < numberRows; i++) {
2512 double value1 = array[i];
2513 double value2 = array2[i];
2514 assert (value1 == value2);
2515 }
2516 }
2517 delete save;
2518 delete [] check;
2519 return returnCode;
2520#else
2521networkBasis_->updateColumn(regionSparse, regionSparse2, -1);
2522return 1;
2523#endif
2524 }
2525#endif
2526}
2527/* Updates one column (FTRAN) from region2
2528 Tries to do FT update
2529 number returned is negative if no room.
2530 Also updates region3
2531 region1 starts as zero and is zero at end */
2532int
2533ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
2534 CoinIndexedVector * regionSparse2,
2535 CoinIndexedVector * regionSparse3,
2536 bool noPermuteRegion3)
2537{
2538#ifdef CLP_DEBUG
2539 regionSparse1->checkClear();
2540#endif
2541 if (!numberRows())
2542 return 0;
2543 int returnCode = 0;
2544#ifndef SLIM_CLP
2545 if (!networkBasis_) {
2546#endif
2547#ifdef CLP_FACTORIZATION_INSTRUMENT
2548 factorization_instrument(-1);
2549#endif
2550 if (coinFactorizationA_) {
2551 coinFactorizationA_->setCollectStatistics(true);
2552 if (coinFactorizationA_->spaceForForrestTomlin()) {
2553 assert (regionSparse2->packedMode());
2554 assert (!regionSparse3->packedMode());
2555 returnCode = coinFactorizationA_->updateTwoColumnsFT(regionSparse1,
2556 regionSparse2,
2557 regionSparse3,
2558 noPermuteRegion3);
2559 } else {
2560 returnCode = coinFactorizationA_->updateColumnFT(regionSparse1,
2561 regionSparse2);
2562 coinFactorizationA_->updateColumn(regionSparse1,
2563 regionSparse3,
2564 noPermuteRegion3);
2565 }
2566 coinFactorizationA_->setCollectStatistics(false);
2567 } else {
2568#if 0
2569 CoinSimpFactorization * fact =
2570 dynamic_cast< CoinSimpFactorization*>(coinFactorizationB_);
2571 if (!fact) {
2572 returnCode = coinFactorizationB_->updateColumnFT(regionSparse1,
2573 regionSparse2);
2574 coinFactorizationB_->updateColumn(regionSparse1,
2575 regionSparse3,
2576 noPermuteRegion3);
2577 } else {
2578 returnCode = fact->updateTwoColumnsFT(regionSparse1,
2579 regionSparse2,
2580 regionSparse3,
2581 noPermuteRegion3);
2582 }
2583#else
2584#ifdef CLP_REUSE_ETAS
2585 int tempInfo[2];
2586 tempInfo[0] = model_->numberIterations();
2587 tempInfo[1] = model_->sequenceIn();
2588 coinFactorizationB_->setUsefulInformation(tempInfo, 3);
2589#endif
2590 returnCode =
2591 coinFactorizationB_->updateTwoColumnsFT(
2592 regionSparse1,
2593 regionSparse2,
2594 regionSparse3,
2595 noPermuteRegion3);
2596#endif
2597 }
2598#ifdef CLP_FACTORIZATION_INSTRUMENT
2599 factorization_instrument(9);
2600#endif
2601#ifdef PRINT_VECTOR
2602 printf("UpdateTwoFT\n");
2603 regionSparse2->print();
2604 regionSparse3->print();
2605#endif
2606 return returnCode;
2607#ifndef SLIM_CLP
2608 } else {
2609 returnCode = updateColumnFT(regionSparse1, regionSparse2);
2610 updateColumn(regionSparse1, regionSparse3, noPermuteRegion3);
2611 }
2612#endif
2613 return returnCode;
2614}
2615/* Updates one column (FTRAN) from region2
2616 number returned is negative if no room
2617 region1 starts as zero and is zero at end */
2618int
2619ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
2620 CoinIndexedVector * regionSparse2,
2621 bool noPermute) const
2622{
2623 if (!noPermute)
2624 regionSparse->checkClear();
2625 if (!coinFactorizationA_->numberRows())
2626 return 0;
2627 coinFactorizationA_->setCollectStatistics(false);
2628 int returnCode = coinFactorizationA_->updateColumn(regionSparse,
2629 regionSparse2,
2630 noPermute);
2631 return returnCode;
2632}
2633/* Updates one column (BTRAN) from region2
2634 region1 starts as zero and is zero at end */
2635int
2636ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
2637 CoinIndexedVector * regionSparse2) const
2638{
2639 if (!numberRows())
2640 return 0;
2641#ifndef SLIM_CLP
2642 if (!networkBasis_) {
2643#endif
2644#ifdef CLP_FACTORIZATION_INSTRUMENT
2645 factorization_instrument(-1);
2646#endif
2647 int returnCode;
2648
2649 if (coinFactorizationA_) {
2650 coinFactorizationA_->setCollectStatistics(true);
2651 returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
2652 regionSparse2);
2653 coinFactorizationA_->setCollectStatistics(false);
2654 } else {
2655 returnCode = coinFactorizationB_->updateColumnTranspose(regionSparse,
2656 regionSparse2);
2657 }
2658#ifdef CLP_FACTORIZATION_INSTRUMENT
2659 factorization_instrument(6);
2660#endif
2661#ifdef PRINT_VECTOR
2662 printf("UpdateTranspose\n");
2663 regionSparse2->print();
2664#endif
2665 return returnCode;
2666#ifndef SLIM_CLP
2667 } else {
2668#ifdef CHECK_NETWORK
2669 CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2670 double * check = new double[coinFactorizationA_->numberRows()];
2671 int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
2672 regionSparse2);
2673 networkBasis_->updateColumnTranspose(regionSparse, save);
2674 int i;
2675 double * array = regionSparse2->denseVector();
2676 int * indices = regionSparse2->getIndices();
2677 int n = regionSparse2->getNumElements();
2678 memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double));
2679 double * array2 = save->denseVector();
2680 int * indices2 = save->getIndices();
2681 int n2 = save->getNumElements();
2682 assert (n == n2);
2683 if (save->packedMode()) {
2684 for (i = 0; i < n; i++) {
2685 check[indices[i]] = array[i];
2686 }
2687 for (i = 0; i < n; i++) {
2688 double value2 = array2[i];
2689 assert (check[indices2[i]] == value2);
2690 }
2691 } else {
2692 int numberRows = coinFactorizationA_->numberRows();
2693 for (i = 0; i < numberRows; i++) {
2694 double value1 = array[i];
2695 double value2 = array2[i];
2696 assert (value1 == value2);
2697 }
2698 }
2699 delete save;
2700 delete [] check;
2701 return returnCode;
2702#else
2703return networkBasis_->updateColumnTranspose(regionSparse, regionSparse2);
2704#endif
2705 }
2706#endif
2707}
2708/* makes a row copy of L for speed and to allow very sparse problems */
2709void
2710ClpFactorization::goSparse()
2711{
2712#ifndef SLIM_CLP
2713 if (!networkBasis_) {
2714#endif
2715 if (coinFactorizationA_) {
2716#ifdef CLP_FACTORIZATION_INSTRUMENT
2717 factorization_instrument(-1);
2718#endif
2719 coinFactorizationA_->goSparse();
2720#ifdef CLP_FACTORIZATION_INSTRUMENT
2721 factorization_instrument(7);
2722#endif
2723 }
2724 }
2725}
2726// Cleans up i.e. gets rid of network basis
2727void
2728ClpFactorization::cleanUp()
2729{
2730#ifndef SLIM_CLP
2731 delete networkBasis_;
2732 networkBasis_ = NULL;
2733#endif
2734 if (coinFactorizationA_)
2735 coinFactorizationA_->resetStatistics();
2736}
2737/// Says whether to redo pivot order
2738bool
2739ClpFactorization::needToReorder() const
2740{
2741#ifdef CHECK_NETWORK
2742 return true;
2743#endif
2744#ifndef SLIM_CLP
2745 if (!networkBasis_)
2746#endif
2747 return true;
2748#ifndef SLIM_CLP
2749 else
2750 return false;
2751#endif
2752}
2753// Get weighted row list
2754void
2755ClpFactorization::getWeights(int * weights) const
2756{
2757#ifdef CLP_FACTORIZATION_INSTRUMENT
2758 factorization_instrument(-1);
2759#endif
2760#ifndef SLIM_CLP
2761 if (networkBasis_) {
2762 // Network - just unit
2763 int numberRows = coinFactorizationA_->numberRows();
2764 for (int i = 0; i < numberRows; i++)
2765 weights[i] = 1;
2766 return;
2767 }
2768#endif
2769 int * numberInRow = coinFactorizationA_->numberInRow();
2770 int * numberInColumn = coinFactorizationA_->numberInColumn();
2771 int * permuteBack = coinFactorizationA_->pivotColumnBack();
2772 int * indexRowU = coinFactorizationA_->indexRowU();
2773 const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
2774 const CoinBigIndex * startRowL = coinFactorizationA_->startRowL();
2775 int numberRows = coinFactorizationA_->numberRows();
2776 if (!startRowL || !coinFactorizationA_->numberInRow()) {
2777 int * temp = new int[numberRows];
2778 memset(temp, 0, numberRows * sizeof(int));
2779 int i;
2780 for (i = 0; i < numberRows; i++) {
2781 // one for pivot
2782 temp[i]++;
2783 CoinBigIndex j;
2784 for (j = startColumnU[i]; j < startColumnU[i] + numberInColumn[i]; j++) {
2785 int iRow = indexRowU[j];
2786 temp[iRow]++;
2787 }
2788 }
2789 CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL();
2790 int * indexRowL = coinFactorizationA_->indexRowL();
2791 int numberL = coinFactorizationA_->numberL();
2792 CoinBigIndex baseL = coinFactorizationA_->baseL();
2793 for (i = baseL; i < baseL + numberL; i++) {
2794 CoinBigIndex j;
2795 for (j = startColumnL[i]; j < startColumnL[i+1]; j++) {
2796 int iRow = indexRowL[j];
2797 temp[iRow]++;
2798 }
2799 }
2800 for (i = 0; i < numberRows; i++) {
2801 int number = temp[i];
2802 int iPermute = permuteBack[i];
2803 weights[iPermute] = number;
2804 }
2805 delete [] temp;
2806 } else {
2807 int i;
2808 for (i = 0; i < numberRows; i++) {
2809 int number = startRowL[i+1] - startRowL[i] + numberInRow[i] + 1;
2810 //number = startRowL[i+1]-startRowL[i]+1;
2811 //number = numberInRow[i]+1;
2812 int iPermute = permuteBack[i];
2813 weights[iPermute] = number;
2814 }
2815 }
2816#ifdef CLP_FACTORIZATION_INSTRUMENT
2817 factorization_instrument(8);
2818#endif
2819}
2820// Set tolerances to safer of existing and given
2821void
2822ClpFactorization::saferTolerances ( double zeroValue,
2823 double pivotValue)
2824{
2825 double newValue;
2826 // better to have small tolerance even if slower
2827 if (zeroValue > 0.0)
2828 newValue = zeroValue;
2829 else
2830 newValue = -zeroTolerance() * zeroValue;
2831 zeroTolerance(CoinMin(zeroTolerance(), zeroValue));
2832 // better to have large tolerance even if slower
2833 if (pivotValue > 0.0)
2834 newValue = pivotValue;
2835 else
2836 newValue = -pivotTolerance() * pivotValue;
2837 pivotTolerance(CoinMin(CoinMax(pivotTolerance(), newValue), 0.999));
2838}
2839// Sets factorization
2840void
2841ClpFactorization::setFactorization(ClpFactorization & rhs)
2842{
2843 ClpFactorization::operator=(rhs);
2844}
2845#endif
2846