1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 #include <assert.h>
17 #include <cstdio>
18 #include <iostream>
19 #include <iomanip>
20 #include <sstream>
21
22 #include "soplex/spxdefines.h"
23 #include "soplex/didxset.h"
24 #include "soplex/mpsinput.h"
25 #include "soplex/spxout.h"
26 #include "soplex/exceptions.h"
27
28 namespace soplex
29 {
30
31 template <class R> typename SPxBasisBase<R>::Desc::Status
32 SPxBasisBase<R>::dualStatus(const SPxColId& id) const
33 {
34 return dualColStatus(static_cast<SPxLPBase<R>*>(theLP)->number(id));
35 }
36
37 template <class R>
38 typename SPxBasisBase<R>::Desc::Status
39 SPxBasisBase<R>::dualStatus(const SPxRowId& id) const
40 {
41 return dualRowStatus((static_cast<SPxLPBase<R>*>(theLP))->number(id));
42 }
43
44 template <class R>
45 typename SPxBasisBase<R>::Desc::Status
46 SPxBasisBase<R>::dualRowStatus(int i) const
47 {
48 assert(theLP != 0);
49
50 if(theLP->rhs(i) < R(infinity))
51 {
52 if(theLP->lhs(i) > R(-infinity))
53 {
54 if(theLP->lhs(i) == theLP->rhs(i))
55 return Desc::D_FREE;
56 else
57 return Desc::D_ON_BOTH;
58 }
59 else
60 return Desc::D_ON_LOWER;
61 }
62 else if(theLP->lhs(i) > R(-infinity))
63 return Desc::D_ON_UPPER;
64 else
65 return Desc::D_UNDEFINED;
66 }
67
68 template <class R>
69 typename SPxBasisBase<R>::Desc::Status
70 SPxBasisBase<R>::dualColStatus(int i) const
71 {
72 assert(theLP != 0);
73
74 if(theLP->SPxLPBase<R>::upper(i) < R(infinity))
75 {
76 if(theLP->SPxLPBase<R>::lower(i) > R(-infinity))
77 {
78 if(theLP->SPxLPBase<R>::lower(i) == theLP->SPxLPBase<R>::upper(i))
79 return Desc::D_FREE;
80 else
81 return Desc::D_ON_BOTH;
82 }
83 else
84 return Desc::D_ON_LOWER;
85 }
86 else if(theLP->SPxLPBase<R>::lower(i) > R(-infinity))
87 return Desc::D_ON_UPPER;
88 else
89 return Desc::D_UNDEFINED;
90 }
91
92 template <class R>
93 void SPxBasisBase<R>::loadMatrixVecs()
94 {
95 assert(theLP != 0);
96 assert(theLP->dim() == matrix.size());
97
98 MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS01 loadMatrixVecs() invalidates factorization"
99 << std::endl;)
100
101 int i;
102 nzCount = 0;
103
104 for(i = theLP->dim() - 1; i >= 0; --i)
105 {
106 matrix[i] = &theLP->vector(baseId(i));
107 nzCount += matrix[i]->size();
108 }
109
110 matrixIsSetup = true;
111 factorized = false;
112
113 if(factor != 0)
114 factor->clear();
115 }
116
117 template <class R>
118 bool SPxBasisBase<R>::isDescValid(const Desc& ds)
119 {
120
121 assert(status() > NO_PROBLEM);
122 assert(theLP != 0);
123
124 int basisdim;
125
126 if(ds.nRows() != theLP->nRows() || ds.nCols() != theLP->nCols())
127 {
128 MSG_DEBUG(std::cout << "IBASIS20 Dimension mismatch\n");
129 return false;
130 }
131
132 basisdim = 0;
133
134 for(int row = ds.nRows() - 1; row >= 0; --row)
135 {
136 if(ds.rowstat[row] >= 0)
137 {
138 if(ds.rowstat[row] != dualRowStatus(row))
139 {
140 MSG_DEBUG(std::cout << "IBASIS21 Basic row " << row << " with incorrect dual status " <<
141 dualRowStatus(row) << "\n");
142 return false;
143 }
144 }
145 else
146 {
147 basisdim++;
148
149 if((ds.rowstat[row] == Desc::P_FIXED
150 && theLP->SPxLPBase<R>::lhs(row) != theLP->SPxLPBase<R>::rhs(row))
151 || (ds.rowstat[row] == Desc::P_ON_UPPER && theLP->SPxLPBase<R>::rhs(row) >= R(infinity))
152 || (ds.rowstat[row] == Desc::P_ON_LOWER && theLP->SPxLPBase<R>::lhs(row) <= R(-infinity)))
153 {
154 MSG_DEBUG(std::cout << "IBASIS22 Nonbasic row with incorrect status: lhs=" <<
155 theLP->SPxLPBase<R>::lhs(row) << ", rhs=" << theLP->SPxLPBase<R>::rhs(
156 row) << ", stat=" << ds.rowstat[row] << "\n");
157 return false;
158 }
159 }
160 }
161
162 for(int col = ds.nCols() - 1; col >= 0; --col)
163 {
164 if(ds.colstat[col] >= 0)
165 {
166 if(ds.colstat[col] != dualColStatus(col))
167 {
168 MSG_DEBUG(std::cout << "IBASIS23 Basic column " << col << " with incorrect dual status " <<
169 ds.colstat[col] << " != " << dualColStatus(col) << "\n");
170 return false;
171 }
172 }
173 else
174 {
175 basisdim++;
176
177 if((ds.colstat[col] == Desc::P_FIXED
178 && theLP->SPxLPBase<R>::lower(col) != theLP->SPxLPBase<R>::upper(col))
179 || (ds.colstat[col] == Desc::P_ON_UPPER && theLP->SPxLPBase<R>::upper(col) >= R(infinity))
180 || (ds.colstat[col] == Desc::P_ON_LOWER && theLP->SPxLPBase<R>::lower(col) <= R(-infinity)))
181 {
182 MSG_DEBUG(std::cout << "IBASIS24 Nonbasic column " << col << " with incorrect status: lower=" <<
183 theLP->SPxLPBase<R>::lower(col) << ", upper=" << theLP->SPxLPBase<R>::upper(
184 col) << ", stat=" << ds.colstat[col] << "\n");
185 return false;
186 }
187 }
188 }
189
190 if(basisdim != theLP->nCols())
191 {
192 MSG_DEBUG(std::cout << "IBASIS25 Incorrect basis dimension " << basisdim << " != " << theLP->nCols()
193 << "\n");
194 return false;
195 }
196
197 // basis descriptor valid
198 return true;
199 }
200
201
202 /*
203 Loading a #Desc# into the basis can be done more efficiently, by
204 explicitely programming both cases, for the rowwise and for the columnwise
205 representation. This implementation hides this distinction in the use of
206 methods #isBasic()# and #vector()#.
207 */
208 template <class R>
209 void SPxBasisBase<R>::loadDesc(const Desc& ds)
210 {
211 assert(status() > NO_PROBLEM);
212 assert(theLP != 0);
213 assert(ds.nRows() == theLP->nRows());
214 assert(ds.nCols() == theLP->nCols());
215
216 SPxId none;
217 int i;
218 int j;
219 bool consistent = true;
220
221 MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS02 loading of Basis invalidates factorization"
222 << std::endl;)
223
224 lastin = none;
225 lastout = none;
226 lastidx = -1;
227 iterCount = 0;
228 updateCount = 0;
229
230 if(&ds != &thedesc)
231 {
232 thedesc = ds;
233 setRep();
234 }
235
236 assert(theLP->dim() == matrix.size());
237
238 nzCount = 0;
239
240 for(j = i = 0; i < theLP->nRows(); ++i)
241 {
242 /* for columns and rows with D_... status, the correct D_... status depends on bounds and sides; if a basis
243 * descriptor is loaded after changing bounds or sides, e.g. in the refine() method, we have to correct them
244 */
245 if(thedesc.rowStatus(i) >= 0)
246 thedesc.rowStatus(i) = dualRowStatus(i);
247 else if(thedesc.rowStatus(i) == SPxBasisBase<R>::Desc::P_FIXED
248 && theLP->SPxLPBase<R>::lhs(i) != theLP->SPxLPBase<R>::rhs(i))
249 {
250 if(theLP->SPxLPBase<R>::lhs(i) > R(-infinity) && theLP->SPxLPBase<R>::maxRowObj(i) < 0.0)
251 thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER;
252 else if(theLP->SPxLPBase<R>::rhs(i) < R(infinity))
253 thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER;
254 else
255 thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_FREE;
256 }
257
258 if(theLP->isBasic(thedesc.rowStatus(i)))
259 {
260 assert(theLP->dim() == matrix.size());
261 assert(j <= matrix.size());
262
263 if(j == matrix.size())
264 {
265 // too many basic variables
266 consistent = false;
267 break;
268 }
269
270 SPxRowId id = theLP->SPxLPBase<R>::rId(i);
271 theBaseId[j] = id;
272 matrix[j] = &theLP->vector(id);
273 nzCount += matrix[j++]->size();
274 }
275 }
276
277 for(i = 0; i < theLP->nCols(); ++i)
278 {
279 /* for columns and rows with D_... status, the correct D_... status depends on bounds and sides; if a basis
280 * descriptor is loaded after changing bounds or sides, e.g. in the refine() method, we have to correct them
281 */
282 if(thedesc.colStatus(i) >= 0)
283 thedesc.colStatus(i) = dualColStatus(i);
284 else if(thedesc.colStatus(i) == SPxBasisBase<R>::Desc::P_FIXED
285 && theLP->SPxLPBase<R>::lower(i) != theLP->SPxLPBase<R>::upper(i))
286 {
287 if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity) && theLP->SPxLPBase<R>::upper(i) >= R(infinity))
288 thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_FREE;
289 else if(theLP->SPxLPBase<R>::upper(i) >= R(infinity)
290 || (theLP->SPxLPBase<R>::lower(i) > R(-infinity) && theLP->SPxLPBase<R>::maxObj(i) < 0.0))
291 thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER;
292 else
293 thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER;
294 }
295
296 if(theLP->isBasic(thedesc.colStatus(i)))
297 {
298 assert(theLP->dim() == matrix.size());
299 assert(j <= matrix.size());
300
301 if(j == matrix.size())
302 {
303 // too many basic variables
304 consistent = false;
305 break;
306 }
307
308 SPxColId id = theLP->SPxLPBase<R>::cId(i);
309 theBaseId[j] = id;
310 matrix[j] = &theLP->vector(id);
311 nzCount += matrix[j++]->size();
312 }
313 }
314
315 if(j < matrix.size())
316 consistent = false; // not enough basic variables
317
318 /* if dimensions are inconsistent, restore slack basis
319 * if dimensions are consistent, then we have setup a correct matrix
320 */
321 if(!consistent)
322 restoreInitialBasis();
323 else
324 matrixIsSetup = true;
325
326 assert(isDescValid(thedesc));
327
328 factorized = false;
329
330 if(factor != 0)
331 factor->clear();
332 }
333
334 template <class R>
335 void SPxBasisBase<R>::setRep()
336 {
337 assert(theLP != 0);
338
339 reDim();
340 minStab = 0.0;
341
342 if(theLP->rep() == SPxSolverBase<R>::ROW)
343 {
344 thedesc.stat = &thedesc.rowstat;
345 thedesc.costat = &thedesc.colstat;
346 }
347 else
348 {
349 thedesc.stat = &thedesc.colstat;
350 thedesc.costat = &thedesc.rowstat;
351 }
352 }
353
354 template <class R>
355 void SPxBasisBase<R>::load(SPxSolverBase<R>* lp, bool initSlackBasis)
356 {
357 assert(lp != 0);
358 theLP = lp;
359
360 setOutstream(*theLP->spxout);
361
362 setRep();
363
364 if(initSlackBasis)
365 {
366 restoreInitialBasis();
367 loadDesc(thedesc);
368 }
369 }
370
371 template <class R>
372 void SPxBasisBase<R>::loadBasisSolver(SLinSolver<R>* p_solver, const bool destroy)
373 {
374 assert(!freeSlinSolver || factor != 0);
375
376 setOutstream(*p_solver->spxout);
377
378 MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS03 loading of Solver invalidates factorization"
379 << std::endl;)
380
381 if(freeSlinSolver)
382 {
383 delete factor;
384 factor = 0;
385 }
386
387 factor = p_solver;
388 factorized = false;
389 factor->clear();
390 freeSlinSolver = destroy;
391 }
392
393 /**
394 * The specification is taken from the
395 *
396 * ILOG CPLEX 7.0 Reference Manual, Appendix E, Page 543.
397 *
398 * This routine should read valid BAS format files.
399 *
400 * @return true if the file was read correctly.
401 *
402 * Here is a very brief outline of the format:
403 *
404 * The format is in a form similar to an MPS file. The basic assumption is that all (column)
405 * variables are nonbasic at their lower bound and all row (variables) are basic; only the
406 * differences to this rule are given. Each data line contains an indicator, a variable name and
407 * possibly a row/constraint name. The following meaning applies with respect to the indicators:
408 *
409 * - XU: the variable is basic, the row is nonbasic at its upper bound
410 * - XL: the variable is basic, the row is nonbasic at its lower bound
411 * - UL: the variable is nonbasic and at its upper bound
412 * - LL: the variable is nonbasic and at its lower bound
413 *
414 * The CPLEX format contains an additional indicator 'BS', but this is unsupported here.
415 *
416 * Nonbasic variables without lower bound have the following default status for SoPlex:
417 * - at their upper bound if finite,
418 * - at zero if free.
419 */
420 template <class R>
421 bool SPxBasisBase<R>::readBasis(
422 std::istream& is,
423 const NameSet* rowNames,
424 const NameSet* colNames)
425 {
426 assert(theLP != 0);
427
428 /* prepare names */
429 const NameSet* rNames = rowNames;
430 const NameSet* cNames = colNames;
431
432 NameSet* p_colNames = 0;
433 NameSet* p_rowNames = 0;
434
435 if(colNames == 0)
436 {
437 int nCols = theLP->nCols();
438 std::stringstream name;
439
440 spx_alloc(p_colNames);
441 p_colNames = new(p_colNames) NameSet();
442 p_colNames->reMax(nCols);
443
444 for(int j = 0; j < nCols; ++j)
445 {
446 name << "x" << j;
447 DataKey key = theLP->colId(j);
448 p_colNames->add(key, name.str().c_str());
449 }
450
451 cNames = p_colNames;
452 }
453
454 if(rNames == 0)
455 {
456 int nRows = theLP->nRows();
457 std::stringstream name;
458
459 spx_alloc(p_rowNames);
460 p_rowNames = new(p_rowNames) NameSet();
461 p_rowNames->reMax(nRows);
462
463 for(int i = 0; i < nRows; ++i)
464 {
465 name << "C" << i;
466 DataKey key = theLP->rowId(i);
467 p_rowNames->add(key, name.str().c_str());
468 }
469
470 rNames = p_rowNames;
471 }
472
473 /* load default basis if necessary */
474 if(status() == NO_PROBLEM)
475 load(theLP, false);
476
477 /* initialize with standard settings */
478 Desc l_desc(thedesc);
479
480 for(int i = 0; i < theLP->nRows(); i++)
481 l_desc.rowstat[i] = dualRowStatus(i);
482
483 for(int i = 0; i < theLP->nCols(); i++)
484 {
485 if(theLP->SPxLPBase<R>::lower(i) == theLP->SPxLPBase<R>::upper(i))
486 l_desc.colstat[i] = Desc::P_FIXED;
487 else if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity)
488 && theLP->SPxLPBase<R>::upper(i) >= R(infinity))
489 l_desc.colstat[i] = Desc::P_FREE;
490 else if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity))
491 l_desc.colstat[i] = Desc::P_ON_UPPER;
492 else
493 l_desc.colstat[i] = Desc::P_ON_LOWER;
494 }
495
496 MPSInput mps(is);
497
498 if(mps.readLine() && (mps.field0() != 0) && !strcmp(mps.field0(), "NAME"))
499 {
500 while(mps.readLine())
501 {
502 int c = -1;
503 int r = -1;
504
505 if((mps.field0() != 0) && !strcmp(mps.field0(), "ENDATA"))
506 {
507 mps.setSection(MPSInput::ENDATA);
508 break;
509 }
510
511 if((mps.field1() == 0) || (mps.field2() == 0))
512 break;
513
514 if((c = cNames->number(mps.field2())) < 0)
515 break;
516
517 if(*mps.field1() == 'X')
518 if(mps.field3() == 0 || (r = rNames->number(mps.field3())) < 0)
519 break;
520
521 if(!strcmp(mps.field1(), "XU"))
522 {
523 l_desc.colstat[c] = dualColStatus(c);
524
525 if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::GREATER_EQUAL)
526 l_desc.rowstat[r] = Desc::P_ON_LOWER;
527 else if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::EQUAL)
528 l_desc.rowstat[r] = Desc::P_FIXED;
529 else
530 l_desc.rowstat[r] = Desc::P_ON_UPPER;
531 }
532 else if(!strcmp(mps.field1(), "XL"))
533 {
534 l_desc.colstat[c] = dualColStatus(c);
535
536 if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::LESS_EQUAL)
537 l_desc.rowstat[r] = Desc::P_ON_UPPER;
538 else if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::EQUAL)
539 l_desc.rowstat[r] = Desc::P_FIXED;
540 else
541 l_desc.rowstat[r] = Desc::P_ON_LOWER;
542 }
543 else if(!strcmp(mps.field1(), "UL"))
544 {
545 l_desc.colstat[c] = Desc::P_ON_UPPER;
546 }
547 else if(!strcmp(mps.field1(), "LL"))
548 {
549 l_desc.colstat[c] = Desc::P_ON_LOWER;
550 }
551 else
552 {
553 mps.syntaxError();
554 break;
555 }
556 }
557 }
558
559 if(!mps.hasError())
560 {
561 if(mps.section() == MPSInput::ENDATA)
562 {
563 // force basis to be different from NO_PROBLEM
564 // otherwise the basis will be overwritten at later stages.
565 setStatus(REGULAR);
566 loadDesc(l_desc);
567 }
568 else
569 mps.syntaxError();
570 }
571
572 if(rowNames == 0)
573 {
574 p_rowNames->~NameSet();
575 spx_free(p_rowNames);
576 }
577
578 if(colNames == 0)
579 {
580 p_colNames->~NameSet();
581 spx_free(p_colNames);
582 }
583
584 #ifndef NDEBUG
585 MSG_DEBUG(thedesc.dump());
586 #endif
587
588 return !mps.hasError();
589 }
590
591
592 /* Get row name - copied from spxmpswrite.cpp
593 *
594 * @todo put this in a common file and unify accross different formats (mps, lp, basis).
595 */
596 template <class R>
597 static const char* getRowName(
598 const SPxLPBase<R>* lp,
599 int idx,
600 const NameSet* rnames,
601 char* buf)
602 {
603 assert(buf != 0);
604 assert(idx >= 0);
605 assert(idx < lp->nRows());
606
607 if(rnames != 0)
608 {
609 DataKey key = lp->rId(idx);
610
611 if(rnames->has(key))
612 return (*rnames)[key];
613 }
614
615 spxSnprintf(buf, 16, "C%d", idx);
616
617 return buf;
618 }
619
620 /* Get column name - copied from spxmpswrite.cpp
621 *
622 * @todo put this in a common file and unify accross different formats (mps, lp, basis).
623 */
624 template <class R>
625 static const char* getColName(
626 const SPxLPBase<R>* lp,
627 int idx,
628 const NameSet* cnames,
629 char* buf)
630 {
631 assert(buf != 0);
632 assert(idx >= 0);
633 assert(idx < lp->nCols());
634
635 if(cnames != 0)
636 {
637 DataKey key = lp->cId(idx);
638
639 if(cnames->has(key))
640 return (*cnames)[key];
641 }
642
643 spxSnprintf(buf, 16, "x%d", idx);
644
645 return buf;
646 }
647
648 /* writes a file in MPS basis format to \p os.
649 *
650 * See SPxBasisBase<R>::readBasis() for a short description of the format.
651 */
652 template <class R>
653 void SPxBasisBase<R>::writeBasis(
654 std::ostream& os,
655 const NameSet* rowNames,
656 const NameSet* colNames,
657 const bool cpxFormat
658 ) const
659 {
660 assert(theLP != 0);
661
662 os.setf(std::ios::left);
663 os << "NAME soplex.bas\n";
664
665 /* do not write basis if there is none */
666 if(status() == NO_PROBLEM)
667 {
668 os << "ENDATA" << std::endl;
669 return;
670 }
671
672 /* start writing */
673 char buf[255];
674 int row = 0;
675
676 for(int col = 0; col < theLP->nCols(); col++)
677 {
678 if(thedesc.colStatus(col) > 0)
679 {
680 /* Find non basic row */
681 for(; row < theLP->nRows(); row++)
682 {
683 if(thedesc.rowStatus(row) < 0)
684 break;
685 }
686
687 assert(row != theLP->nRows());
688
689 if(thedesc.rowStatus(row) == Desc::P_ON_UPPER && (!cpxFormat
690 || theLP->LPRowSetBase<R>::type(row) == LPRowBase<R>::RANGE))
691 os << " XU ";
692 else
693 os << " XL ";
694
695 os << std::setw(8) << getColName(theLP, col, colNames, buf);
696
697 /* break in two parts since buf is reused */
698 os << " "
699 << getRowName(theLP, row, rowNames, buf)
700 << std::endl;
701
702 row++;
703 }
704 else
705 {
706 if(thedesc.colStatus(col) == Desc::P_ON_UPPER)
707 {
708 os << " UL "
709 << getColName(theLP, col, colNames, buf)
710 << std::endl;
711 }
712 else
713 {
714 /* Default is all non-basic variables on lower bound (if finite) or at zero (if free).
715 * nothing to do in this case.
716 */
717 assert(theLP->lower(col) <= R(-infinity) || thedesc.colStatus(col) == Desc::P_ON_LOWER
718 || thedesc.colStatus(col) == Desc::P_FIXED);
719 assert(theLP->lower(col) > R(-infinity) || theLP->upper(col) < R(infinity)
720 || thedesc.colStatus(col) == Desc::P_FREE);
721 }
722 }
723 }
724
725 #ifndef NDEBUG
726 MSG_DEBUG(thedesc.dump());
727
728 // Check that we covered all nonbasic rows - the remaining should be basic.
729 for(; row < theLP->nRows(); row++)
730 {
731 if(thedesc.rowStatus(row) < 0)
732 break;
733 }
734
735 assert(row == theLP->nRows());
736
737 #endif // NDEBUG
738
739 os << "ENDATA" << std::endl;
740 }
741
742 template <class R>
743 void SPxBasisBase<R>::printMatrix() const
744 {
745
746 assert(matrixIsSetup);
747
748 for(int i = 0; i < matrix.size(); i++)
749 {
750 std::cout << "C" << i << "=" << *matrix[i] << std::endl;
751 }
752 }
753
754 template <class R>
755 void SPxBasisBase<R>::printMatrixMTX(int number)
756 {
757 int dim;
758 int nnz;
759 char filename[SPX_MAXSTRLEN];
760
761 dim = matrix.size();
762 nnz = nzCount;
763 spxSnprintf(filename, SPX_MAXSTRLEN, "basis/basis%d.mtx", number);
764 std::cout << "printing basis matrix to file " << filename << "\n";
765 FILE* basisfile;
766 basisfile = fopen(filename, "w");
767 // print marker necessary for reading the file in Matlab
768 fprintf(basisfile, "%%%%MatrixMarket matrix coordinate Real general\n");
769 // print matrix information
770 fprintf(basisfile, "%d %d %d\n", dim, dim, nnz);
771
772 // print matrix data
773 for(int i = 0; i < matrix.size(); ++i)
774 {
775 for(int j = 0; j < baseVec(i).size(); ++j)
776 {
777 int idx = baseVec(i).index(j);
778 R val = baseVec(i).value(j);
779 fprintf(basisfile, "%d %d %.13" REAL_FORMAT "\n", i + 1, idx + 1, val);
780 }
781 }
782
783 fclose(basisfile);
784
785 return;
786 }
787
788 template <class R>
789 void SPxBasisBase<R>::change(
790 int i,
791 SPxId& id,
792 const SVectorBase<R>* enterVec,
793 const SSVectorBase<R>* eta)
794 {
795
796 assert(matrixIsSetup);
797 assert(!id.isValid() || (enterVec != 0));
798 assert(factor != 0);
799
800 lastidx = i;
801 lastin = id;
802
803 if(id.isValid() && i >= 0)
804 {
805 assert(enterVec != 0);
806
807 // update the counter for nonzeros in the basis matrix
808 nzCount = nzCount - matrix[i]->size() + enterVec->size();
809 // let the new id enter the basis
810 matrix[i] = enterVec;
811 lastout = theBaseId[i];
812 theBaseId[i] = id;
813
814 ++iterCount;
815 ++updateCount;
816
817 MSG_DEBUG(std::cout << "factor_stats: iteration= " << this->iteration()
818 << " update= " << updateCount
819 << " total_update= " << totalUpdateCount
820 << " nonzero_B= " << nzCount
821 << " nonzero_LU= " << factor->memory()
822 << " factor_fill= " << lastFill
823 << " time= " << theLP->time()
824 << std::endl;)
825
826 // never factorize? Just do it !
827 if(!factorized)
828 factorize();
829
830 // too much memory growth ?
831 else if(R(factor->memory()) > 1000 + factor->dim() + lastMem * memFactor)
832 {
833 MSG_INFO3((*this->spxout), (*this->spxout) <<
834 "IBASIS04 memory growth factor triggers refactorization"
835 << " memory= " << factor->memory()
836 << " lastMem= " << lastMem
837 << " memFactor= " << memFactor
838 << std::endl;)
839 factorize();
840 }
841
842 // relative fill too high ?
843 else if(R(factor->memory()) > lastFill * R(nzCount))
844 {
845 MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS04 fill factor triggers refactorization"
846 << " memory= " << factor->memory()
847 << " nzCount= " << nzCount
848 << " lastFill= " << lastFill
849 << std::endl;)
850
851 factorize();
852 }
853 // absolute fill in basis matrix too high ?
854 else if(nzCount > lastNzCount)
855 {
856 MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS05 nonzero factor triggers refactorization"
857 << " nzCount= " << nzCount
858 << " lastNzCount= " << lastNzCount
859 << " nonzeroFactor= " << nonzeroFactor
860 << std::endl;)
861 factorize();
862 }
863 // too many updates ?
864 else if(updateCount >= maxUpdates)
865 {
866 MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS06 update count triggers refactorization"
867 << " updateCount= " << updateCount
868 << " maxUpdates= " << maxUpdates
869 << std::endl;)
870 factorize();
871 }
872 else
873 {
874 try
875 {
876 #ifdef MEASUREUPDATETIME
877 theTime.start();
878 #endif
879 factor->change(i, *enterVec, eta);
880 totalUpdateCount++;
881 #ifdef MEASUREUPDATETIME
882 theTime.stop();
883 #endif
884 }
885 catch(...)
886 {
887 MSG_INFO3((*this->spxout), (*this->spxout) <<
888 "IBASIS13 problems updating factorization; refactorizing basis"
889 << std::endl;)
890
891 #ifdef MEASUREUPDATETIME
892 theTime.stop();
893 #endif
894
895 // singularity was detected in update; we refactorize
896 factorize();
897
898 // if factorize() detects singularity, an exception is thrown, hence at this point we have a regular basis
899 // and can try the update again
900 assert(status() >= SPxBasisBase<R>::REGULAR);
901
902 try
903 {
904 #ifdef MEASUREUPDATETIME
905 theTime.start();
906 #endif
907 factor->change(i, *enterVec, eta);
908 totalUpdateCount++;
909 #ifdef MEASUREUPDATETIME
910 theTime.stop();
911 #endif
912 }
913 // with a freshly factorized, regular basis, the update is unlikely to fail; if this happens nevertheless,
914 // we have to invalidate the basis to have the statuses correct
915 catch(const SPxException& F)
916 {
917 MSG_INFO3((*this->spxout), (*this->spxout) <<
918 "IBASIS14 problems updating factorization; invalidating factorization"
919 << std::endl;)
920
921 #ifdef MEASUREUPDATETIME
922 theTime.stop();
923 #endif
924
925 factorized = false;
926 throw F;
927 }
928 }
929
930 assert(minStab > 0.0);
931
932 if(factor->status() != SLinSolver<R>::OK || factor->stability() < minStab)
933 {
934 MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS07 stability triggers refactorization"
935 << " stability= " << factor->stability()
936 << " minStab= " << minStab
937 << std::endl;)
938 factorize();
939 }
940 }
941 }
942 else
943 lastout = id;
944 }
945
946 template <class R>
947 void SPxBasisBase<R>::factorize()
948 {
949
950 assert(factor != 0);
951
952 if(!matrixIsSetup)
953 loadDesc(thedesc);
954
955 assert(matrixIsSetup);
956
957 updateCount = 0;
958
959 switch(factor->load(matrix.get_ptr(), matrix.size()))
960 {
961 case SLinSolver<R>::OK :
962 if(status() == SINGULAR)
963 setStatus(REGULAR);
964
965 factorized = true;
966 minStab = factor->stability();
967
968 // This seems always to be about 1e-7
969 if(minStab > 1e-4)
970 minStab *= 0.001;
971
972 if(minStab > 1e-5)
973 minStab *= 0.01;
974
975 if(minStab > 1e-6)
976 minStab *= 0.1;
977
978 break;
979
980 case SLinSolver<R>::SINGULAR :
981 setStatus(SINGULAR);
982 factorized = false;
983 break;
984
985 default :
986 MSG_ERROR(std::cerr << "EBASIS08 error: unknown status of factorization.\n";)
987 factorized = false;
988 throw SPxInternalCodeException("XBASIS01 This should never happen.");
989 }
990
991 // get nonzero count of factorization
992 lastMem = factor->memory();
993 // compute fill ratio between factorization and basis matrix (multiplied with tolerance)
994 lastFill = fillFactor * R(lastMem) / R(nzCount > 0 ? nzCount : 1);
995 lastNzCount = int(nonzeroFactor * R(nzCount > 0 ? nzCount : 1));
996
997 if(status() == SINGULAR)
998 {
999 throw SPxStatusException("Cannot factorize singular matrix");
1000 }
1001 }
1002
1003 template <class R>
1004 VectorBase<R>& SPxBasisBase<R>::multWithBase(VectorBase<R>& x) const
1005 {
1006 assert(status() > SINGULAR);
1007 assert(theLP->dim() == x.dim());
1008
1009 int i;
1010 VectorBase<R> tmp(x);
1011
1012 if(!matrixIsSetup)
1013 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1014
1015 assert(matrixIsSetup);
1016
1017 for(i = x.dim() - 1; i >= 0; --i)
1018 x[i] = *(matrix[i]) * tmp;
1019
1020 return x;
1021 }
1022
1023 template <class R>
1024 void SPxBasisBase<R>::multWithBase(SSVectorBase<R>& x, SSVectorBase<R>& result) const
1025 {
1026 assert(status() > SINGULAR);
1027 assert(theLP->dim() == x.dim());
1028 assert(x.dim() == result.dim());
1029
1030 if(!matrixIsSetup)
1031 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1032
1033 result.clear();
1034
1035 assert(matrixIsSetup);
1036
1037 for(int i = 0; i < x.dim(); ++i)
1038 result.add(i, (*matrix[i]) * x);
1039
1040 return;
1041 }
1042
1043 template <class R>
1044 VectorBase<R>& SPxBasisBase<R>::multBaseWith(VectorBase<R>& x) const
1045 {
1046 assert(status() > SINGULAR);
1047 assert(theLP->dim() == x.dim());
1048
1049 int i;
1050 VectorBase<R> tmp(x);
1051
1052 if(!matrixIsSetup)
1053 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1054
1055 assert(matrixIsSetup);
1056
1057 x.clear();
1058
1059 for(i = x.dim() - 1; i >= 0; --i)
1060 {
1061 if(tmp[i] != 0.0)
1062 x.multAdd(tmp[i], *(matrix[i]));
1063 }
1064
1065 return x;
1066 }
1067
1068 template <class R>
1069 void SPxBasisBase<R>::multBaseWith(SSVectorBase<R>& x, SSVectorBase<R>& result) const
1070 {
1071 assert(status() > SINGULAR);
1072 assert(theLP->dim() == x.dim());
1073 assert(x.dim() == result.dim());
1074
1075 if(!matrixIsSetup)
1076 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1077
1078 assert(matrixIsSetup);
1079
1080 result.clear();
1081
1082 if(x.isSetup())
1083 {
1084 for(int i = 0; i < x.size(); ++i)
1085 {
1086 int idx = x.index(i);
1087 result.multAdd(x[idx], (*matrix[idx]));
1088 }
1089 }
1090 else
1091 {
1092 for(int i = 0; i < x.dim(); ++i)
1093 result.multAdd(x[i], (*matrix[i]));
1094 }
1095
1096 return;
1097 }
1098
1099 template <class R>
1100 /* compute an estimated condition number for the current basis matrix
1101 * by computing estimates of the norms of B and B^-1 using the power method.
1102 * maxiters and tolerance control the accuracy of the estimate.
1103 */
1104 R SPxBasisBase<R>::condition(int maxiters, R tolerance)
1105 {
1106 int dimension = matrix.size();
1107 int miniters = 3; // minimum number of power method iterations
1108 int i;
1109 int c;
1110 R norm;
1111 R norminv;
1112 R norm1;
1113 R norm2;
1114
1115 // catch corner case of empty matrix
|
(1) Event cond_false: |
Condition "dimension == 0", taking false branch. |
1116 if(dimension == 0)
|
(2) Event if_end: |
End of if statement. |
1117 return 1.0;
1118
|
(3) Event write_zero_model: |
"SSVectorBase" sets "x.idx" to "nullptr". [details] |
| Also see events: |
[var_deref_model] |
1119 SSVectorBase<R> x(dimension);
1120 SSVectorBase<R> y(dimension);
1121
1122 // check whether a regular basis matrix is available
|
(4) Event cond_false: |
Condition "this->status() < soplex::SPxBasisBase<double>::REGULAR", taking false branch. |
1123 if(status() < REGULAR)
|
(5) Event if_end: |
End of if statement. |
1124 return 0;
1125
|
(6) Event cond_false: |
Condition "!this->matrixIsSetup", taking false branch. |
1126 if(!matrixIsSetup)
|
(7) Event if_end: |
End of if statement. |
1127 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1128
|
(8) Event cond_false: |
Condition "!this->factorized", taking false branch. |
1129 if(!factorized)
|
(9) Event if_end: |
End of if statement. |
1130 factorize();
1131
1132 // initialize vectors
1133 norm1 = 1.0 / (R) dimension;
1134
|
(10) Event cond_true: |
Condition "i < dimension", taking true branch. |
1135 for(i = 0; i < dimension; i++)
1136 x.add(i, norm1);
1137
1138 y = x;
1139
1140 // compute norm of B
1141 for(c = 0; c < maxiters; ++c)
1142 {
1143 norm2 = norm1;
1144
1145 // y = B*x
1146 multBaseWith(x, y);
1147 norm1 = y.length();
1148
1149 // stop if converged
1150 if(c >= miniters && spxAbs(norm1 - norm2) < tolerance * norm1)
1151 break;
1152
1153 // x = B^T*y and normalize
1154 multWithBase(y, x);
1155 norm2 = 1.0 / x.length();
1156 x *= norm2;
1157 }
1158
1159 norm = norm1;
1160
1161 // reinitialize vectors
1162 x.clear();
1163 y.clear();
1164 norm1 = 1.0 / (R) dimension;
1165
1166 for(i = 0; i < dimension; i++)
1167 x.add(i, norm1);
1168
1169 y = x;
1170
1171 // compute norm of B^-1
1172 for(c = 0; c < maxiters; ++c)
1173 {
1174 norm2 = norm1;
1175
1176 // x = B^-1*y
1177 factor->solveRight(x, y);
1178 x.setup();
1179 norm1 = x.length();
1180
1181 // stop if converged
1182 if(c >= miniters && spxAbs(norm1 - norm2) < tolerance * norm1)
1183 break;
1184
1185 // y = B^-T*x and normalize
1186 factor->solveLeft(y, x);
1187 y.setup();
1188 norm2 = 1.0 / y.length();
1189 y *= norm2;
1190 }
1191
1192 norminv = norm1;
1193
1194 return norm * norminv;
1195 }
1196
1197 /* compute one of several matrix metrics based on the diagonal of the LU factorization */
1198 template <class R>
1199 R SPxBasisBase<R>::getMatrixMetric(int type)
1200 {
1201 R metric = R(infinity);
1202
1203 if(factorized)
1204 metric = factor->matrixMetric(type);
1205
1206 return metric;
1207 }
1208
1209 template <class R>
1210 void SPxBasisBase<R>::dump()
1211 {
1212 assert(status() > NO_PROBLEM);
1213 assert(theLP != 0);
1214 assert(thedesc.nRows() == theLP->nRows());
1215 assert(thedesc.nCols() == theLP->nCols());
1216 assert(theLP->dim() == matrix.size());
1217
1218 int i, basesize;
1219
1220 // Dump regardless of the verbosity level if this method is called.
1221
1222 std::cout << "DBASIS09 Basis entries:";
1223 basesize = 0;
1224
1225 for(i = 0; i < theLP->nRows(); ++i)
1226 {
1227 if(theLP->isBasic(thedesc.rowStatus(i)))
1228 {
1229 if(basesize % 10 == 0)
1230 std::cout << std::endl << "DBASIS10 ";
1231
1232 SPxRowId id = theLP->SPxLPBase<R>::rId(i);
1233 std::cout << "\tR" << theLP->number(id);
1234 basesize++;
1235 }
1236 }
1237
1238 for(i = 0; i < theLP->nCols(); ++i)
1239 {
1240 if(theLP->isBasic(thedesc.colStatus(i)))
1241 {
1242 if(basesize % 10 == 0)
1243 std::cout << std::endl << "DBASIS11 ";
1244
1245 SPxColId id = theLP->SPxLPBase<R>::cId(i);
1246 std::cout << "\tC" << theLP->number(id);
1247 basesize++;
1248 }
1249 }
1250
1251 std::cout << std::endl;
1252
1253 assert(basesize == matrix.size());
1254 }
1255
1256 template <class R>
1257
1258 bool SPxBasisBase<R>::isConsistent() const
1259 {
1260 #ifdef ENABLE_CONSISTENCY_CHECKS
1261 int primals = 0;
1262 int i;
1263
1264 if(status() > NO_PROBLEM)
1265 {
1266 if(theLP == 0)
1267 return MSGinconsistent("SPxBasisBase<R>");
1268
1269 if(theBaseId.size() != theLP->dim() || matrix.size() != theLP->dim())
1270 return MSGinconsistent("SPxBasisBase<R>");
1271
1272 if(thedesc.nCols() != theLP->nCols() || thedesc.nRows() != theLP->nRows())
1273 return MSGinconsistent("SPxBasisBase<R>");
1274
1275 for(i = 0; i < thedesc.nRows(); ++i)
1276 {
1277 if(thedesc.rowStatus(i) >= 0)
1278 {
1279 if(thedesc.rowStatus(i) != dualRowStatus(i))
1280 return MSGinconsistent("SPxBasisBase<R>");
1281 }
1282 else
1283 ++primals;
1284 }
1285
1286 for(i = 0; i < thedesc.nCols(); ++i)
1287 {
1288 if(thedesc.colStatus(i) >= 0)
1289 {
1290 if(thedesc.colStatus(i) != dualColStatus(i))
1291 return MSGinconsistent("SPxBasisBase<R>");
1292 }
1293 else
1294 ++primals;
1295 }
1296
1297 if(primals != thedesc.nCols())
1298 return MSGinconsistent("SPxBasisBase<R>");
1299 }
1300
1301 return thedesc.isConsistent() && theBaseId.isConsistent()
1302 && matrix.isConsistent() && factor->isConsistent();
1303 #else
1304 return true;
1305 #endif // CONSISTENCY_CHECKS
1306 }
1307
1308 template <class R>
1309 SPxBasisBase<R>::SPxBasisBase(Timer::TYPE ttype)
1310 : theLP(0)
1311 , matrixIsSetup(false)
1312 , factor(0)
1313 , factorized(false)
1314 , maxUpdates(200)
1315 , nonzeroFactor(10.0)
1316 , fillFactor(5.0)
1317 , memFactor(1.5)
1318 , iterCount(0)
1319 , lastIterCount(0)
1320 , iterDegenCheck(0)
1321 , updateCount(0)
1322 , totalUpdateCount(0)
1323 , nzCount(1)
1324 , lastMem(0)
1325 , lastFill(0)
1326 , lastNzCount(0)
1327 , theTime(0)
1328 , timerType(ttype)
1329 , lastidx(0)
1330 , minStab(0.0)
1331 , thestatus(NO_PROBLEM)
1332 , freeSlinSolver(false)
1333 , spxout(0)
1334 {
1335 // info: is not consistent at this moment, e.g. because theLP == 0
1336
1337 theTime = TimerFactory::createTimer(timerType);
1338 }
1339
1340
1341 /**@warning Do not change the LP object.
1342 * Only pointer to that object is copied.
1343 * Hint: no problem, we use this function for copy
1344 * constructor of SPxSolverBase<R>
1345 */
1346 template <class R>
1347 SPxBasisBase<R>::SPxBasisBase(const SPxBasisBase<R>& old)
1348 : theLP(old.theLP)
1349 , theBaseId(old.theBaseId)
1350 , matrix(old.matrix)
1351 , matrixIsSetup(old.matrixIsSetup)
1352 , factor(old.factor)
1353 , factorized(old.factorized)
1354 , maxUpdates(old.maxUpdates)
1355 , nonzeroFactor(old.nonzeroFactor)
1356 , fillFactor(old.fillFactor)
1357 , memFactor(old.memFactor)
1358 , iterCount(old.iterCount)
1359 , lastIterCount(old.lastIterCount)
1360 , iterDegenCheck(old.iterDegenCheck)
1361 , updateCount(old.updateCount)
1362 , totalUpdateCount(old.totalUpdateCount)
1363 , nzCount(old.nzCount)
1364 , lastMem(old.lastMem)
1365 , lastFill(old.lastFill)
1366 , lastNzCount(old.lastNzCount)
1367 , theTime(old.theTime)
1368 , lastin(old.lastin)
1369 , lastout(old.lastout)
1370 , lastidx(old.lastidx)
1371 , minStab(old.minStab)
1372 , thestatus(old.thestatus)
1373 , thedesc(old.thedesc)
1374 , spxout(old.spxout)
1375 {
1376 theTime = TimerFactory::createTimer(old.theTime->type());
1377
1378 this->factor = old.factor->clone();
1379 freeSlinSolver = true;
1380
1381 assert(SPxBasisBase<R>::isConsistent());
1382 }
1383
1384 template <class R>
1385 SPxBasisBase<R>::~SPxBasisBase<R>()
1386 {
1387
1388 assert(!freeSlinSolver || factor != 0);
1389
1390 if(freeSlinSolver)
1391 {
1392 delete factor;
1393 factor = 0;
1394 }
1395
1396 theTime->~Timer();
1397 spx_free(theTime);
1398 }
1399
1400 template <class R>
1401
1402 /**@warning Note that we do not create a deep copy of the corresponding SPxSolverBase<R> object.
1403 * Only the reference to that object is copied.
1404 */
1405 SPxBasisBase<R>& SPxBasisBase<R>::operator=(const SPxBasisBase<R>& rhs)
1406 {
1407
1408 assert(!freeSlinSolver || factor != 0);
1409
1410 if(this != &rhs)
1411 {
1412 theLP = rhs.theLP;
1413 theBaseId = rhs.theBaseId;
1414 matrix = rhs.matrix;
1415 matrixIsSetup = rhs.matrixIsSetup;
1416
1417 if(freeSlinSolver)
1418 {
1419 delete factor;
1420 factor = 0;
1421 }
1422
1423 factor = rhs.factor->clone();
1424 freeSlinSolver = true;
1425
1426 factorized = rhs.factorized;
1427 maxUpdates = rhs.maxUpdates;
1428 nonzeroFactor = rhs.nonzeroFactor;
1429 fillFactor = rhs.fillFactor;
1430 memFactor = rhs.memFactor;
1431 iterCount = rhs.iterCount;
1432 nzCount = rhs.nzCount;
1433 lastFill = rhs.lastFill;
1434 lastNzCount = rhs.lastNzCount;
1435 lastin = rhs.lastin;
1436 lastout = rhs.lastout;
1437 lastidx = rhs.lastidx;
1438 minStab = rhs.minStab;
1439 thestatus = rhs.thestatus;
1440 thedesc = rhs.thedesc;
1441
1442 assert(SPxBasisBase<R>::isConsistent());
1443 }
1444
1445 return *this;
1446 }
1447
1448
1449
1450 //
1451 // Auxiliary functions.
1452 //
1453
1454 // Pretty-printing of basis status.
1455 template <class R> // Why can't we remove the class R and make it empty?
1456 std::ostream& operator<<(std::ostream& os,
1457 const typename SPxBasisBase<R>::SPxStatus& status)
1458 {
1459 switch(status)
1460 {
1461 case SPxBasisBase<R>::NO_PROBLEM:
1462 os << "NO_PROBLEM";
1463 break;
1464
1465 case SPxBasisBase<R>::SINGULAR:
1466 os << "SINGULAR";
1467 break;
1468
1469 case SPxBasisBase<R>::REGULAR:
1470 os << "REGULAR";
1471 break;
1472
1473 case SPxBasisBase<R>::DUAL:
1474 os << "DUAL";
1475 break;
1476
1477 case SPxBasisBase<R>::PRIMAL:
1478 os << "PRIMAL";
1479 break;
1480
1481 case SPxBasisBase<R>::OPTIMAL:
1482 os << "OPTIMAL";
1483 break;
1484
1485 case SPxBasisBase<R>::UNBOUNDED:
1486 os << "UNBOUNDED";
1487 break;
1488
1489 case SPxBasisBase<R>::INFEASIBLE:
1490 os << "INFEASIBLE";
1491 break;
1492
1493 default:
1494 os << "?other?";
1495 break;
1496 }
1497
1498 return os;
1499 }
1500
1501
1502 } // namespace soplex
1503