1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2022 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
17 #include <assert.h>
18 #include <iostream>
19
20 #include "soplex/spxdefines.h"
|
(1) Event include_recursion: |
#include file "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scipoptsuite-8.0.1/soplex/src/soplex/spxweightst.h" includes itself: spxweightst.h -> spxweightst.hpp -> spxweightst.h |
|
(2) Event caretline: |
^ |
21 #include "soplex/spxweightst.h"
22 #include "soplex/svset.h"
23 #include "soplex/sorter.h"
24
25 namespace soplex
26 {
27 #define EPS 1e-6
28 #define STABLE 1e-3 // the sparsest row/column may only have a pivot of size STABLE*maxEntry
29
30
31 template <class R>
32 bool SPxWeightST<R>::isConsistent() const
33 {
34 #ifdef ENABLE_CONSISTENCY_CHECKS
35 return rowWeight.isConsistent()
36 && colWeight.isConsistent()
37 && rowRight.isConsistent()
38 && colUp.isConsistent();
39 // && SPxStarter::isConsistent(); // not yet implemented
40 #else
41 return true;
42 #endif
43 }
44
45 /* Generating the Starting Basis
46 The generation of a starting basis works as follows: After setting up the
47 preference arrays #weight# and #coWeight#, #Id#s are selected to be dual in
48 a greedy manner. Initially, the first #Id# is taken. Then the next #Id# is
49 checked wheter its vector is linearly dependend of the vectors of the #Id#s
50 allready selected. If not, it is added to them. This is iterated until a
51 full matrix has been constructed.
52
53 Testing for linear independence is done very much along the lines of LU
54 factorization. A vector is taken, and updated with all previous L-vectors.
55 Then the maximal absolut element is selected as pivot element for computing
56 the L-vector. This is stored for further processing.
57 */
58
59 /*
60 The following two functions set the status of |id| to primal or dual,
61 respectively.
62 */
63 template <class R>
64 void SPxWeightST<R>::setPrimalStatus(
65 typename SPxBasisBase<R>::Desc& desc,
66 const SPxSolverBase<R>& base,
67 const SPxId& id)
68 {
69 if(id.isSPxRowId())
70 {
71 int n = base.number(SPxRowId(id));
72
73 if(base.rhs(n) >= R(infinity))
74 {
75 if(base.lhs(n) <= R(-infinity))
76 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_FREE;
77 else
78 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER;
79 }
80 else
81 {
82 if(base.lhs(n) <= R(-infinity))
83 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER;
84 else if(base.lhs(n) >= base.rhs(n) - base.epsilon())
85 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_FIXED;
86 else if(rowRight[n])
87 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER;
88 else
89 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER;
90 }
91 }
92 else
93 {
94 int n = base.number(SPxColId(id));
95
96 if(base.SPxLPBase<R>::upper(n) >= R(infinity))
97 {
98 if(base.SPxLPBase<R>::lower(n) <= R(-infinity))
99 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_FREE;
100 else
101 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER;
102 }
103 else
104 {
105 if(base.SPxLPBase<R>::lower(n) <= R(-infinity))
106 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER;
107 else if(base.SPxLPBase<R>::lower(n) >= base.SPxLPBase<R>::upper(n) - base.epsilon())
108 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_FIXED;
109 else if(colUp[n])
110 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER;
111 else
112 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER;
113 }
114 }
115 }
116
117 // ----------------------------------------------------------------
118 template <class R>
119 static void setDualStatus(
120 typename SPxBasisBase<R>::Desc& desc,
121 const SPxSolverBase<R>& base,
122 const SPxId& id)
123 {
124 if(id.isSPxRowId())
125 {
126 int n = base.number(SPxRowId(id));
127 desc.rowStatus(n) = base.basis().dualRowStatus(n);
128 }
129 else
130 {
131 int n = base.number(SPxColId(id));
132 desc.colStatus(n) = base.basis().dualColStatus(n);
133 }
134 }
135 // ----------------------------------------------------------------
136
137 /// Compare class for row weights, used for sorting.
138 template <typename T>
139 struct Compare
140 {
141 public:
142 /// constructor
143 Compare() : weight(0) {}
144 // const SPxSolverBase* base; ///< the solver
145 const T* weight; ///< the weights to compare
146
147 /// compares the weights
148 T operator()(int i1, int i2) const
149 {
150 return weight[i1] - weight[i2];
151 }
152 };
153
154 // ----------------------------------------------------------------
155 /**
156 The following method initializes \p pref such that it contains the
157 set of \ref soplex::SPxId "SPxIds" ordered following \p rowWeight and
158 \p colWeight. For the sorting we take the following approach: first
159 we sort the rows, then the columns. Finally we perform a mergesort
160 of both.
161 */
162 template <class R>
163 static void initPrefs(
164 DataArray < SPxId >& pref,
165 const SPxSolverBase<R>& base,
166 const Array<R>& rowWeight,
167 const Array<R>& colWeight)
168 {
169 DataArray<int> row(base.nRows());
170 DataArray<int> col(base.nCols());
171 int i;
172 int j;
173 int k;
174
175 Compare<R> compare;
176 // compare.base = &base;
177
178 for(i = 0; i < base.nRows(); ++i)
179 row[i] = i;
180
181 compare.weight = rowWeight.get_const_ptr();
182
183 SPxQuicksort(row.get_ptr(), row.size(), compare); // Sort rows
184
185 for(i = 0; i < base.nCols(); ++i)
186 col[i] = i;
187
188 compare.weight = colWeight.get_const_ptr();
189
190 SPxQuicksort(col.get_ptr(), col.size(), compare); // Sort column
191
192 i = 0;
193 j = 0;
194 k = 0;
195
196 while(k < pref.size()) // merge sort
197 {
198 if(rowWeight[row[i]] < colWeight[col[j]])
199 {
200 pref[k++] = base.rId(row[i++]);
201
202 if(i >= base.nRows())
203 while(k < pref.size())
204 pref[k++] = base.cId(col[j++]);
205 }
206 else
207 {
208 pref[k++] = base.cId(col[j++]);
209
210 if(j >= base.nCols())
211 while(k < pref.size())
212 pref[k++] = base.rId(row[i++]);
213 }
214 }
215
216 assert(i == base.nRows());
217 assert(j == base.nCols());
218 }
219
220 // ----------------------------------------------------------------
221 template <class R>
222 void SPxWeightST<R>::generate(SPxSolverBase<R>& base)
223 {
224 SPxId tmpId;
225
226 forbidden.reSize(base.dim());
227 rowWeight.reSize(base.nRows());
228 colWeight.reSize(base.nCols());
229 rowRight.reSize(base.nRows());
230 colUp.reSize(base.nCols());
231
232 if(base.rep() == SPxSolverBase<R>::COLUMN)
233 {
234 weight = &colWeight;
235 coWeight = &rowWeight;
236 }
237 else
238 {
239 weight = &rowWeight;
240 coWeight = &colWeight;
241 }
242
243 assert(weight->size() == base.coDim());
244 assert(coWeight->size() == base.dim());
245
246 setupWeights(base);
247
248 typename SPxBasisBase<R>::Desc desc(base);
249 // desc.reSize(base.nRows(), base.nCols());
250
251 DataArray < SPxId > pref(base.nRows() + base.nCols());
252 initPrefs(pref, base, rowWeight, colWeight);
253
254 int i;
255 int stepi;
256 int j;
257 int sel;
258
259 for(i = 0; i < base.dim(); ++i)
260 forbidden[i] = 0;
261
262 if(base.rep() == SPxSolverBase<R>::COLUMN)
263 {
264 // in COLUMN rep we scan from beginning to end
265 i = 0;
266 stepi = 1;
267 }
268 else
269 {
270 // in ROW rep we scan from end to beginning
271 i = pref.size() - 1;
272 stepi = -1;
273 }
274
275 int dim = base.dim();
276 R maxEntry = 0;
277
278 for(; i >= 0 && i < pref.size(); i += stepi)
279 {
280 tmpId = pref[i];
281 const SVectorBase<R>& vec = base.vector(tmpId);
282 sel = -1;
283
284 // column or row singleton ?
285 if(vec.size() == 1)
286 {
287 int idx = vec.index(0);
288
289 if(forbidden[idx] < 2)
290 {
291 sel = idx;
292 dim += (forbidden[idx] > 0) ? 1 : 0;
293 }
294 }
295 else
296 {
297 maxEntry = vec.maxAbs();
298
299 // initialize the nonzero counter
300 int minRowEntries = base.nRows();
301
302 // find a stable index with a sparse row/column
303 for(j = vec.size(); --j >= 0;)
304 {
305 R x = vec.value(j);
306 int k = vec.index(j);
307 int nRowEntries = base.coVector(k).size();
308
309 if(!forbidden[k] && (spxAbs(x) > STABLE * maxEntry) && (nRowEntries < minRowEntries))
310 {
311 minRowEntries = nRowEntries;
312 sel = k;
313 }
314 }
315 }
316
317 // we found a valid index
318 if(sel >= 0)
319 {
320 MSG_DEBUG(
321
322 if(pref[i].type() == SPxId::ROW_ID)
323 std::cout << "DWEIST01 r" << base.number(pref[i]);
324 else
325 std::cout << "DWEIST02 c" << base.number(pref[i]);
326 )
327
328 forbidden[sel] = 2;
329
330 // put current column/row into basis
331 if(base.rep() == SPxSolverBase<R>::COLUMN)
332 setDualStatus(desc, base, pref[i]);
333 else
334 setPrimalStatus(desc, base, pref[i]);
335
336 for(j = vec.size(); --j >= 0;)
337 {
338 R x = vec.value(j);
339 int k = vec.index(j);
340
341 if(!forbidden[k] && (x > EPS * maxEntry || -x > EPS * maxEntry))
342 {
343 forbidden[k] = 1;
344 --dim;
345 }
346 }
347
348 if(--dim == 0)
349 {
350 //@ for(++i; i < pref.size(); ++i)
351 if(base.rep() == SPxSolverBase<R>::COLUMN)
352 {
353 // set all remaining indeces to nonbasic status
354 for(i += stepi; i >= 0 && i < pref.size(); i += stepi)
355 setPrimalStatus(desc, base, pref[i]);
356
357 // fill up the basis wherever linear independence is assured
358 for(i = forbidden.size(); --i >= 0;)
359 {
360 if(forbidden[i] < 2)
361 setDualStatus(desc, base, base.coId(i));
362 }
363 }
364 else
365 {
366 for(i += stepi; i >= 0 && i < pref.size(); i += stepi)
367 setDualStatus(desc, base, pref[i]);
368
369 for(i = forbidden.size(); --i >= 0;)
370 {
371 if(forbidden[i] < 2)
372 setPrimalStatus(desc, base, base.coId(i));
373 }
374 }
375
376 break;
377 }
378 }
379 // sel == -1
380 else if(base.rep() == SPxSolverBase<R>::COLUMN)
381 setPrimalStatus(desc, base, pref[i]);
382 else
383 setDualStatus(desc, base, pref[i]);
384
385 #ifndef NDEBUG
386 {
387 int n, m;
388
389 for(n = 0, m = forbidden.size(); n < forbidden.size(); ++n)
390 m -= (forbidden[n] != 0) ? 1 : 0;
391
392 assert(m == dim);
393 }
394 #endif // NDEBUG
395 }
396
397 assert(dim == 0);
398
399 base.loadBasis(desc);
400 #ifdef TEST
401 base.init();
402
403 int changed = 0;
404 const VectorBase<R>& pvec = base.pVec();
405
406 for(i = pvec.dim() - 1; i >= 0; --i)
407 {
408 if(desc.colStatus(i) == SPxBasisBase<R>::Desc::P_ON_UPPER
409 && base.lower(i) > R(-infinity) && pvec[i] > base.maxObj(i))
410 {
411 changed = 1;
412 desc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER;
413 }
414 else if(desc.colStatus(i) == SPxBasisBase<R>::Desc::P_ON_LOWER
415 && base.upper(i) < R(infinity) && pvec[i] < base.maxObj(i))
416 {
417 changed = 1;
418 desc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER;
419 }
420 }
421
422 if(changed)
423 {
424 std::cout << "changed basis\n";
425 base.loadBasis(desc);
426 }
427 else
428 std::cout << "nothing changed\n";
429
430 #endif // TEST
431 }
432
433 // ----------------------------------------------------------------
434
435 /* Computation of Weights
436 */
437 template <class R>
438 void SPxWeightST<R>::setupWeights(SPxSolverBase<R>& base)
439 {
440 const VectorBase<R>& obj = base.maxObj();
441 const VectorBase<R>& low = base.lower();
442 const VectorBase<R>& up = base.upper();
443 const VectorBase<R>& rhs = base.rhs();
444 const VectorBase<R>& lhs = base.lhs();
445 int i;
446
447 R eps = base.epsilon();
448 R maxabs = 1.0;
449
450 // find absolut biggest entry in bounds and left-/right hand side
451 for(i = 0; i < base.nCols(); i++)
452 {
453 if((up[i] < R(infinity)) && (spxAbs(up[i]) > maxabs))
454 maxabs = spxAbs(up[i]);
455
456 if((low[i] > R(-infinity)) && (spxAbs(low[i]) > maxabs))
457 maxabs = spxAbs(low[i]);
458 }
459
460 for(i = 0; i < base.nRows(); i++)
461 {
462 if((rhs[i] < R(infinity)) && (spxAbs(rhs[i]) > maxabs))
463 maxabs = spxAbs(rhs[i]);
464
465 if((lhs[i] > R(-infinity)) && (spxAbs(lhs[i]) > maxabs))
466 maxabs = spxAbs(lhs[i]);
467 }
468
469 /**@todo The comments are wrong. The first is for dual simplex and
470 * the secound for primal one. Is anything else wrong?
471 * Also the values are nearly the same for both cases.
472 * Should this be ? Changed the values for
473 * r_fixed to 0 because of maros-r7. It is not clear why
474 * this makes a difference because all constraints in that
475 * instance are of equality type.
476 * Why is rowRight sometimes not set?
477 */
478 if(base.rep() * base.type() > 0) // primal simplex
479 {
480 const R ax = 1e-3 / obj.maxAbs();
481 const R bx = 1.0 / maxabs;
482 const R nne = ax / base.nRows(); // 1e-4 * ax;
483 const R c_fixed = 1e+5;
484 const R r_fixed = 0; // TK20010103: was 1e+4 (maros-r7)
485 const R c_dbl_bounded = 1e+1;
486 const R r_dbl_bounded = 0;
487 const R c_bounded = 1e+1;
488 const R r_bounded = 0;
489 const R c_free = -1e+4;
490 const R r_free = -1e+5;
491
492 for(i = base.nCols() - 1; i >= 0; i--)
493 {
494 R n = nne * (base.colVector(i).size() - 1); // very small value that is zero for col singletons
495 R x = ax * obj[i]; // this is at most 1e-3, probably a lot smaller
496 R u = bx * up [i]; // this is at most 1, probably a lot smaller
497 R l = bx * low[i]; // this is at most 1, probably a lot smaller
498
499 if(up[i] < R(infinity))
500 {
501 if(spxAbs(low[i] - up[i]) < eps)
502 colWeight[i] = c_fixed + n + spxAbs(x);
503 else if(low[i] > R(-infinity))
504 {
505 colWeight[i] = c_dbl_bounded + l - u + n;
506
507 l = spxAbs(l);
508 u = spxAbs(u);
509
510 if(u < l)
511 {
512 colUp[i] = true;
513 colWeight[i] += x;
514 }
515 else
516 {
517 colUp[i] = false;
518 colWeight[i] -= x;
519 }
520 }
521 else
522 {
523 colWeight[i] = c_bounded - u + x + n;
524 colUp[i] = true;
525 }
526 }
527 else
528 {
529 if(low[i] > R(-infinity))
530 {
531 colWeight[i] = c_bounded + l + n - x;
532 colUp[i] = false;
533 }
534 else
535 {
536 colWeight[i] = c_free + n - spxAbs(x);
537 }
538 }
539 }
540
541 for(i = base.nRows() - 1; i >= 0; i--)
542 {
543 if(rhs[i] < R(infinity))
544 {
545 if(spxAbs(lhs[i] - rhs[i]) < eps)
546 {
547 rowWeight[i] = r_fixed;
548 }
549 else if(lhs[i] > R(-infinity))
550 {
551 R u = bx * rhs[i];
552 R l = bx * lhs[i];
553
554 rowWeight[i] = r_dbl_bounded + l - u;
555 rowRight[i] = spxAbs(u) < spxAbs(l);
556 }
557 else
558 {
559 rowWeight[i] = r_bounded - bx * rhs[i];
560 rowRight[i] = true;
561 }
562 }
563 else
564 {
565 if(lhs[i] > R(-infinity))
566 {
567 rowWeight[i] = r_bounded + bx * lhs[i];
568 rowRight[i] = false;
569 }
570 else
571 {
572 rowWeight[i] = r_free;
573 }
574 }
575 }
576 }
577 else
578 {
579 assert(base.rep() * base.type() < 0); // dual simplex
580
581 const R ax = 1.0 / obj.maxAbs();
582 const R bx = 1e-2 / maxabs;
583 const R nne = 1e-4 * bx;
584 const R c_fixed = 1e+5;
585 const R r_fixed = 1e+4;
586 const R c_dbl_bounded = 1;
587 const R r_dbl_bounded = 0;
588 const R c_bounded = 0;
589 const R r_bounded = 0;
590 const R c_free = -1e+4;
591 const R r_free = -1e+5;
592
593 for(i = base.nCols() - 1; i >= 0; i--)
594 {
595 R n = nne * (base.colVector(i).size() - 1);
596 R x = ax * obj[i];
597 R u = bx * up [i];
598 R l = bx * low[i];
599
600 if(up[i] < R(infinity))
601 {
602 if(spxAbs(low[i] - up[i]) < eps)
603 colWeight[i] = c_fixed + n + spxAbs(x);
604 else if(low[i] > R(-infinity))
605 {
606 if(x > 0)
607 {
608 colWeight[i] = c_dbl_bounded + x - u + n;
609 colUp[i] = true;
610 }
611 else
612 {
613 colWeight[i] = c_dbl_bounded - x + l + n;
614 colUp[i] = false;
615 }
616 }
617 else
618 {
619 colWeight[i] = c_bounded - u + x + n;
620 colUp[i] = true;
621 }
622 }
623 else
624 {
625 if(low[i] > R(-infinity))
626 {
627 colWeight[i] = c_bounded - x + l + n;
628 colUp[i] = false;
629 }
630 else
631 colWeight[i] = c_free + n - spxAbs(x);
632 }
633 }
634
635 for(i = base.nRows() - 1; i >= 0; i--)
636 {
637 const R len1 = 1; // (base.rowVector(i).length() + base.epsilon());
638 R n = 0; // nne * (base.rowVector(i).size() - 1);
639 R u = bx * len1 * rhs[i];
640 R l = bx * len1 * lhs[i];
641 R x = ax * len1 * (obj * base.rowVector(i));
642
643 if(rhs[i] < R(infinity))
644 {
645 if(spxAbs(lhs[i] - rhs[i]) < eps)
646 rowWeight[i] = r_fixed + n + spxAbs(x);
647 else if(lhs[i] > R(-infinity))
648 {
649 if(x > 0)
650 {
651 rowWeight[i] = r_dbl_bounded + x - u + n;
652 rowRight[i] = true;
653 }
654 else
655 {
656 rowWeight[i] = r_dbl_bounded - x + l + n;
657 rowRight[i] = false;
658 }
659 }
660 else
661 {
662 rowWeight[i] = r_bounded - u + n + x;
663 rowRight[i] = true;
664 }
665 }
666 else
667 {
668 if(lhs[i] > R(-infinity))
669 {
670 rowWeight[i] = r_bounded + l + n - x;
671 rowRight[i] = false;
672 }
673 else
674 {
675 rowWeight[i] = r_free + n - spxAbs(x);
676 }
677 }
678 }
679 }
680
681 MSG_DEBUG(
682 {
683 for(i = 0; i < base.nCols(); i++)
684 std::cout << "C i= " << i
685 << " up= " << colUp[i]
686 << " w= " << colWeight[i]
687 << std::endl;
688 for(i = 0; i < base.nRows(); i++)
689 std::cout << "R i= " << i
690 << " rr= " << rowRight[i]
691 << " w= " << rowWeight[i]
692 << std::endl;
693 })
694 }
695 } // namespace soplex
696