1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24
25 #include <assert.h>
26 #include "soplex/cring.h"
27
28 namespace soplex
29 {
30 /* This number is used to decide wether a value is zero
31 * or was explicitly set to zero.
32 */
33 /*
34 #define SOPLEX_MARKER 1e-100
35 */
36
37 static const Real verySparseFactorRat = 0.001;
38 static const Real verySparseFactor4rightRat = 0.2;
39
40 /* generic heap management */
41 static void enQueueMaxRat(int* heap, int* size, int elem)
42 {
43 int i, j;
44
45 j = (*size)++;
46
47 while(j > 0)
48 {
49 i = (j - 1) / 2;
50
51 if(elem > heap[i])
52 {
53 heap[j] = heap[i];
54 j = i;
55 }
56 else
57 break;
58 }
59
60 heap[j] = elem;
61
62 #ifdef SOPLEX_DEBUG
63
64 // no NDEBUG define, since this block is really expensive
65 for(i = 1; i < *size; ++i)
66 for(j = 0; j < i; ++j)
67 assert(heap[i] != heap[j]);
68
69 #endif /* SOPLEX_DEBUG */
70 }
71
72 static int deQueueMaxRat(int* heap, int* size)
73 {
74 int e, elem;
75 int i, j, s;
76 int e1, e2;
77
78 elem = *heap;
79 e = heap[s = --(*size)];
80 --s;
81
82 for(j = 0, i = 1; i < s; i = 2 * j + 1)
83 {
84 e1 = heap[i];
85 e2 = heap[i + 1];
86
87 if(e1 > e2)
88 {
89 if(e < e1)
90 {
91 heap[j] = e1;
92 j = i;
93 }
94 else
95 {
96 heap[j] = e;
97 return elem;
98 }
99 }
100 else
101 {
102 if(e < e2)
103 {
104 heap[j] = e2;
105 j = i + 1;
106 }
107 else
108 {
109 heap[j] = e;
110 return elem;
111 }
112 }
113 }
114
115 if(i < *size && e < heap[i])
116 {
117 heap[j] = heap[i];
118 j = i;
119 }
120
121 heap[j] = e;
122
123 return elem;
124 }
125
126 static void enQueueMinRat(int* heap, int* size, int elem)
127 {
128 int i, j;
129
130 j = (*size)++;
131
132 while(j > 0)
133 {
134 i = (j - 1) / 2;
135
136 if(elem < heap[i])
137 {
138 heap[j] = heap[i];
139 j = i;
140 }
141 else
142 break;
143 }
144
145 heap[j] = elem;
146
147 #ifdef SOPLEX_DEBUG
148
149 // no NDEBUG define, since this block is really expensive
150 for(i = 1; i < *size; ++i)
151 for(j = 0; j < i; ++j)
152 assert(heap[i] != heap[j]);
153
154 #endif /* SOPLEX_DEBUG */
155 }
156
157 static int deQueueMinRat(int* heap, int* size)
158 {
159 int e, elem;
160 int i, j, s;
161 int e1, e2;
162
163 elem = *heap;
164 e = heap[s = --(*size)];
165 --s;
166
167 for(j = 0, i = 1; i < s; i = 2 * j + 1)
168 {
169 e1 = heap[i];
170 e2 = heap[i + 1];
171
172 if(e1 < e2)
173 {
174 if(e > e1)
175 {
176 heap[j] = e1;
177 j = i;
178 }
179 else
180 {
181 heap[j] = e;
182 return elem;
183 }
184 }
185 else
186 {
187 if(e > e2)
188 {
189 heap[j] = e2;
190 j = i + 1;
191 }
192 else
193 {
194 heap[j] = e;
195 return elem;
196 }
197 }
198 }
199
200 if(i < *size && e > heap[i])
201 {
202 heap[j] = heap[i];
203 j = i;
204 }
205
206 heap[j] = e;
207
208 return elem;
209 }
210
211 /************************************************************/
212 inline CLUFactorRational::Temp::Temp()
213 : s_mark(0)
214 , s_cact(0)
215 , stage(0)
216 , pivot_col(0)
217 , pivot_colNZ(0)
218 , pivot_row(0)
219 , pivot_rowNZ(0)
220 {}
221
222 inline void CLUFactorRational::Temp::init(int p_dim)
223 {
224 s_max.reDim(p_dim);
225 spx_realloc(s_cact, p_dim);
226 spx_realloc(s_mark, p_dim);
227 stage = 0;
228 }
229
230 inline void CLUFactorRational::Temp::clear()
231 {
232 if(s_mark != 0)
233 spx_free(s_mark);
234
235 if(s_cact != 0)
236 spx_free(s_cact);
237
238 if(pivot_col != 0)
239 spx_free(pivot_col);
240
241 if(pivot_colNZ != 0)
242 spx_free(pivot_colNZ);
243
244 if(pivot_row != 0)
245 spx_free(pivot_row);
246
247 if(pivot_rowNZ != 0)
248 spx_free(pivot_rowNZ);
249
250 try
251 {
252 s_max.reDim(0);
253 }
254 catch(const SPxMemoryException& x)
255 {
(1) Event exception_thrown: |
An exception of type "soplex::SPxMemoryException const" is thrown. |
256 throw x;
257 }
258 }
259
(1) Event exn_spec_violation: |
An exception of type "soplex::SPxMemoryException const" is thrown but the exception specification "/*implicit*/noexcept" doesn't allow it to be thrown. This will result in a call to terminate(). |
Also see events: |
[fun_call_w_exception] |
260 inline CLUFactorRational::Temp::~Temp()
261 {
(2) Event fun_call_w_exception: |
Called function throws an exception of type "soplex::SPxMemoryException const". [details] |
Also see events: |
[exn_spec_violation] |
262 clear();
263 }
264
265 /************************************************************/
266 inline void CLUFactorRational::initPerm()
267 {
268
269 for(int i = 0; i < thedim; ++i)
270 row.orig[i] = row.perm[i] = col.orig[i] = col.perm[i] = -1;
271 }
272
273 /*****************************************************************************/
274
275 inline void CLUFactorRational::setPivot(const int p_stage,
276 const int p_col,
277 const int p_row,
278 const Rational& val)
279 {
280 assert(row.perm[p_row] < 0);
281 assert(col.perm[p_col] < 0);
282
283 row.orig[p_stage] = p_row;
284 col.orig[p_stage] = p_col;
285 row.perm[p_row] = p_stage;
286 col.perm[p_col] = p_stage;
287 diag[p_row] = 1.0 / val;
288
289 if(spxAbs(diag[p_row]) > maxabs)
290 maxabs = spxAbs(diag[p_row]);
291 }
292
293 /*****************************************************************************/
294 /*
295 * Perform garbage collection on row file
296 */
297 inline void CLUFactorRational::packRows()
298 {
299 int n, i, j, l_row;
300 Dring* ring, *list;
301
302 int* l_ridx = u.row.idx;
303 VectorRational& l_rval = u.row.val;
304 int* l_rlen = u.row.len;
305 int* l_rmax = u.row.max;
306 int* l_rbeg = u.row.start;
307
308 n = 0;
309 list = &(u.row.list);
310
311 for(ring = list->next; ring != list; ring = ring->next)
312 {
313 l_row = ring->idx;
314
315 if(l_rbeg[l_row] != n)
316 {
317 do
318 {
319 l_row = ring->idx;
320 i = l_rbeg[l_row];
321 assert(l_rlen[l_row] <= l_rmax[l_row]);
322 l_rbeg[l_row] = n;
323 l_rmax[l_row] = l_rlen[l_row];
324 j = i + l_rlen[l_row];
325
326 for(; i < j; ++i, ++n)
327 {
328 assert(n <= i);
329 l_ridx[n] = l_ridx[i];
330 l_rval[n] = l_rval[i];
331 }
332
333 ring = ring->next;
334 }
335 while(ring != list);
336
337 goto terminatePackRows;
338 }
339
340 n += l_rlen[l_row];
341
342 l_rmax[l_row] = l_rlen[l_row];
343 }
344
345 terminatePackRows:
346
347 u.row.max[thedim] = 0;
348 u.row.used = n;
349 }
350
351 /*****************************************************************************/
352 /*
353 * Perform garbage collection on column file
354 */
355 inline void CLUFactorRational::forestPackColumns()
356 {
357 int n, i, j, colno;
358 Dring* ring, *list;
359
360 VectorRational& cval = u.col.val;
361 int* cidx = u.col.idx;
362 int* clen = u.col.len;
363 int* cmax = u.col.max;
364 int* cbeg = u.col.start;
365
366 n = 0;
367 list = &u.col.list;
368
369 for(ring = list->next; ring != list; ring = ring->next)
370 {
371 colno = ring->idx;
372
373 if(cbeg[colno] != n)
374 {
375 do
376 {
377 colno = ring->idx;
378 i = cbeg[colno];
379 cbeg[colno] = n;
380 cmax[colno] = clen[colno];
381 j = i + clen[colno];
382
383 for(; i < j; ++i)
384 {
385 cval[n] = cval[i];
386 cidx[n++] = cidx[i];
387 }
388
389 ring = ring->next;
390 }
391 while(ring != list);
392
393 goto terminatePackColumns;
394 }
395
396 n += clen[colno];
397
398 cmax[colno] = clen[colno];
399 }
400
401 terminatePackColumns :
402
403 u.col.used = n;
404 u.col.max[thedim] = 0;
405 }
406
407 /*
408 * Make row of fac large enough to hold len nonzeros.
409 */
410 inline void CLUFactorRational::remaxRow(int p_row, int len)
411 {
412 assert(u.row.max[p_row] < len);
413
414 if(u.row.elem[p_row].next == &(u.row.list)) /* last in row file */
415 {
416 int delta = len - u.row.max[p_row];
417
418 if(delta > u.row.val.dim() - u.row.used)
419 {
420 packRows();
421 delta = len - u.row.max[p_row]; // packRows() changes u.row.max[] !
422
423 if(u.row.val.dim() < rowMemMult * u.row.used + len)
424 minRowMem(2 * u.row.used + len);
425
426 /* minRowMem(rowMemMult * u.row.used + len); */
427 }
428
429 assert(delta <= u.row.val.dim() - u.row.used
430
431 && "ERROR: could not allocate memory for row file");
432
433 u.row.used += delta;
434 u.row.max[p_row] = len;
435 }
436 else /* row must be moved to end of row file */
437 {
438 int i, j, k;
439 VectorRational& val = u.row.val;
440 Dring* ring;
441
442 if(len > u.row.val.dim() - u.row.used)
443 {
444 packRows();
445
446 if(u.row.val.dim() < rowMemMult * u.row.used + len)
447 minRowMem(2 * u.row.used + len);
448
449 /* minRowMem(rowMemMult * u.row.used + len);*/
450 }
451
452 assert(len <= u.row.val.dim() - u.row.used
453
454 && "ERROR: could not allocate memory for row file");
455
456 int* idx = u.row.idx;
457 j = u.row.used;
458 i = u.row.start[p_row];
459 k = u.row.len[p_row] + i;
460 u.row.start[p_row] = j;
461 u.row.used += len;
462
463 u.row.max[u.row.elem[p_row].prev->idx] += u.row.max[p_row];
464 u.row.max[p_row] = len;
465 removeDR(u.row.elem[p_row]);
466 ring = u.row.list.prev;
467 init2DR(u.row.elem[p_row], *ring);
468
469 for(; i < k; ++i, ++j)
470 {
471 val[j] = val[i];
472 idx[j] = idx[i];
473 }
474 }
475
476 assert(u.row.start[u.row.list.prev->idx] + u.row.max[u.row.list.prev->idx]
477
478 == u.row.used);
479 }
480
481 /*************************************************************************/
482 /*
483 * Perform garbage collection on column file
484 */
485 inline void CLUFactorRational::packColumns()
486 {
487 int n, i, j, l_col;
488 Dring* ring, *list;
489
490 int* l_cidx = u.col.idx;
491 int* l_clen = u.col.len;
492 int* l_cmax = u.col.max;
493 int* l_cbeg = u.col.start;
494
495 n = 0;
496 list = &(u.col.list);
497
498 for(ring = list->next; ring != list; ring = ring->next)
499 {
500 l_col = ring->idx;
501
502 if(l_cbeg[l_col] != n)
503 {
504 do
505 {
506 l_col = ring->idx;
507 i = l_cbeg[l_col];
508 l_cbeg[l_col] = n;
509 l_cmax[l_col] = l_clen[l_col];
510 j = i + l_clen[l_col];
511
512 for(; i < j; ++i)
513 l_cidx[n++] = l_cidx[i];
514
515 ring = ring->next;
516 }
517 while(ring != list);
518
519 goto terminatePackColumns;
520 }
521
522 n += l_clen[l_col];
523
524 l_cmax[l_col] = l_clen[l_col];
525 }
526
527 terminatePackColumns :
528
529 u.col.used = n;
530 u.col.max[thedim] = 0;
531 }
532
533 /*
534 * Make column col of fac large enough to hold len nonzeros.
535 */
536 inline void CLUFactorRational::remaxCol(int p_col, int len)
537 {
538 assert(u.col.max[p_col] < len);
539
540 if(u.col.elem[p_col].next == &(u.col.list)) /* last in column file */
541 {
542 int delta = len - u.col.max[p_col];
543
544 if(delta > u.col.size - u.col.used)
545 {
546 packColumns();
547 delta = len - u.col.max[p_col];
548
549 if(u.col.size < colMemMult * u.col.used + len)
550 minColMem(2 * u.col.used + len);
551
552 /* minColMem(colMemMult * u.col.used + len); */
553 }
554
555 assert(delta <= u.col.size - u.col.used
556
557 && "ERROR: could not allocate memory for column file");
558
559 u.col.used += delta;
560 u.col.max[p_col] = len;
561 }
562
563 else /* column must be moved to end of column file */
564 {
565 int i, j, k;
566 int* idx;
567 Dring* ring;
568
569 if(len > u.col.size - u.col.used)
570 {
571 packColumns();
572
573 if(u.col.size < colMemMult * u.col.used + len)
574 minColMem(2 * u.col.used + len);
575
576 /* minColMem(colMemMult * u.col.used + len); */
577 }
578
579 assert(len <= u.col.size - u.col.used
580
581 && "ERROR: could not allocate memory for column file");
582
583 j = u.col.used;
584 i = u.col.start[p_col];
585 k = u.col.len[p_col] + i;
586 u.col.start[p_col] = j;
587 u.col.used += len;
588
589 u.col.max[u.col.elem[p_col].prev->idx] += u.col.max[p_col];
590 u.col.max[p_col] = len;
591 removeDR(u.col.elem[p_col]);
592 ring = u.col.list.prev;
593 init2DR(u.col.elem[p_col], *ring);
594
595 idx = u.col.idx;
596
597 for(; i < k; ++i)
598 idx[j++] = idx[i];
599 }
600 }
601
602 /*
603 * Make column col of fac large enough to hold len nonzeros.
604 */
605 inline void CLUFactorRational::forestReMaxCol(int p_col, int len)
606 {
607 assert(u.col.max[p_col] < len);
608
609 if(u.col.elem[p_col].next == &(u.col.list)) /* last in column file */
610 {
611 int delta = len - u.col.max[p_col];
612
613 if(delta > u.col.size - u.col.used)
614 {
615 forestPackColumns();
616 delta = len - u.col.max[p_col];
617
618 if(u.col.size < colMemMult * u.col.used + len)
619 forestMinColMem(int(colMemMult * u.col.used + len));
620 }
621
622 assert(delta <= u.col.size - u.col.used
623
624 && "ERROR: could not allocate memory for column file");
625
626 u.col.used += delta;
627 u.col.max[p_col] = len;
628 }
629
630 else /* column must be moved to end of column file */
631 {
632 int i, j, k;
633 int* idx = u.col.idx;
634 VectorRational& val = u.col.val;
635 Dring* ring;
636
637 if(len > u.col.size - u.col.used)
638 {
639 forestPackColumns();
640
641 if(u.col.size < colMemMult * u.col.used + len)
642 forestMinColMem(int(colMemMult * u.col.used + len));
643 }
644
645 assert(len <= u.col.size - u.col.used
646
647 && "ERROR: could not allocate memory for column file");
648
649 j = u.col.used;
650 i = u.col.start[p_col];
651 k = u.col.len[p_col] + i;
652 u.col.start[p_col] = j;
653 u.col.used += len;
654
655 u.col.max[u.col.elem[p_col].prev->idx] += u.col.max[p_col];
656 u.col.max[p_col] = len;
657 removeDR(u.col.elem[p_col]);
658 ring = u.col.list.prev;
659 init2DR(u.col.elem[p_col], *ring);
660
661 for(; i < k; ++i)
662 {
663 val[j] = val[i];
664 idx[j++] = idx[i];
665 }
666 }
667 }
668
669 /*****************************************************************************/
670
671 /**
672 * \brief Performs the Forrest-Tomlin update of the LU factorization.
673 *
674 * BH: I suppose this is implemented as described in UH Suhl, LM Suhl: A fast LU
675 * update for linear programming, Annals of OR 43, p. 33-47, 1993.
676 *
677 * @param p_col Index of basis column to replace.
678 * @param p_work Dense vector to substitute in the basis.
679 * @param num Number of nonzeros in vector represented by p_work.
680 * @param nonz Indices of nonzero elements in vector p_work.
681 *
682 * The parameters num and nonz are used for the following optimization: If both
683 * are nonzero, indices of the nonzero elements provided in nonz (num giving
684 * their number) allow to access only those nonzero entries. Otherwise we have
685 * to go through the entire dense vector element by element.
686 *
687 * After copying p_work into U, p_work is used to expand the row r, which is
688 * needed to restore the triangular structure of U.
689 *
690 * Also num and nonz are used to maintain a heap if there are only very few
691 * nonzeros to be eliminated. This is plainly wrong if the method is called with
692 * nonz==0, see todo at the corresponding place below.
693 *
694 * @throw SPxStatusException if the loaded matrix is singular
695 *
696 * @todo Use an extra member variable as a buffer for working with the dense
697 * row instead of misusing p_work. I think that should be as efficient and
698 * much cleaner.
699 */
700
701 inline void CLUFactorRational::forestUpdate(int p_col, Rational* p_work, int num, int* nonz)
702 {
703 int i, j, k, h, m, n;
704 int ll, c, r, rowno;
705 Rational x;
706
707 int* lbeg = l.start;
708
709 VectorRational& cval = u.col.val;
710 int* cidx = u.col.idx;
711 int* cmax = u.col.max;
712 int* clen = u.col.len;
713 int* cbeg = u.col.start;
714
715 VectorRational& rval = u.row.val;
716 int* ridx = u.row.idx;
717 int* rmax = u.row.max;
718 int* rlen = u.row.len;
719 int* rbeg = u.row.start;
720
721 int* rperm = row.perm;
722 int* rorig = row.orig;
723 int* cperm = col.perm;
724 int* corig = col.orig;
725
726 Rational l_maxabs = maxabs;
727 int dim = thedim;
728
729 /* Remove column p_col from U
730 */
731 j = cbeg[p_col];
732 i = clen[p_col];
733 nzCnt -= i;
734
735 for(i += j - 1; i >= j; --i)
736 {
737 m = cidx[i]; // remove column p_col from row m
738 k = rbeg[m];
739 h = --(rlen[m]) + k; // decrease length of row m
740
741 while(ridx[k] != p_col)
742 ++k;
743
744 assert(k <= h); // k is the position of p_col, h is last position
745
746 ridx[k] = ridx[h]; // store last index at the position of p_col
747
748 rval[k] = rval[h];
749 }
750
751 /* Insert new vector column p_col thereby determining the highest permuted
752 * row index r.
753 *
754 * Distinguish between optimized call (num > 0, nonz != 0) and
755 * non-optimized one.
756 */
757 assert(num); // otherwise the assert( nonz != 0 ) below should fail
758
759 if(num)
760 {
761 // Optimized call.
762 assert(nonz != 0);
763
764 clen[p_col] = 0;
765
766 if(num > cmax[p_col])
767 forestReMaxCol(p_col, num);
768
769 cidx = u.col.idx;
770
771 cval = u.col.val;
772
773 k = cbeg[p_col];
774
775 r = 0;
776
777 for(j = 0; j < num; ++j)
778 {
779 i = nonz[j];
780 x = p_work[i];
781 p_work[i] = 0;
782
783 if(x != 0)
784 {
785 if(spxAbs(x) > l_maxabs)
786 l_maxabs = spxAbs(x);
787
788 /* insert to column file */
789 assert(k - cbeg[p_col] < cmax[p_col]);
790
791 cval[k] = x;
792
793 cidx[k++] = i;
794
795 /* insert to row file */
796 if(rmax[i] <= rlen[i])
797 {
798 remaxRow(i, rlen[i] + 1);
799 rval = u.row.val;
800 ridx = u.row.idx;
801 }
802
803 h = rbeg[i] + (rlen[i])++;
804
805 rval[h] = x;
806 ridx[h] = p_col;
807
808 /* check permuted row index */
809
810 if(rperm[i] > r)
811 r = rperm[i];
812 }
813 }
814
815 nzCnt += (clen[p_col] = k - cbeg[p_col]);
816 }
817 else
818 {
819 // Non-optimized call: We have to access all elements of p_work.
820 assert(nonz == 0);
821
822 /*
823 * clen[col] = 0;
824 * reMaxCol(fac, col, dim);
825 */
826 cidx = u.col.idx;
827 cval = u.col.val;
828 k = cbeg[p_col];
829 j = k + cmax[p_col];
830 r = 0;
831
832 for(i = 0; i < dim; ++i)
833 {
834 x = p_work[i];
835 p_work[i] = 0;
836
837 if(x != 0)
838 {
839 if(spxAbs(x) > l_maxabs)
840 l_maxabs = spxAbs(x);
841
842 /* insert to column file */
843 if(k >= j)
844 {
845 clen[p_col] = k - cbeg[p_col];
846 forestReMaxCol(p_col, dim - i);
847 cidx = u.col.idx;
848 cval = u.col.val;
849 k = cbeg[p_col];
850 j = k + cmax[p_col];
851 k += clen[p_col];
852 }
853
854 assert(k - cbeg[p_col] < cmax[p_col]);
855
856 cval[k] = x;
857 cidx[k++] = i;
858
859 /* insert to row file */
860
861 if(rmax[i] <= rlen[i])
862 {
863 remaxRow(i, rlen[i] + 1);
864 rval = u.row.val;
865 ridx = u.row.idx;
866 }
867
868 h = rbeg[i] + (rlen[i])++;
869
870 rval[h] = x;
871 ridx[h] = p_col;
872
873 /* check permuted row index */
874
875 if(rperm[i] > r)
876 r = rperm[i];
877 }
878 }
879
880 nzCnt += (clen[p_col] = k - cbeg[p_col]);
881
882 if(cbeg[p_col] + cmax[p_col] == u.col.used)
883 {
884 u.col.used -= cmax[p_col];
885 cmax[p_col] = clen[p_col];
886 u.col.used += cmax[p_col];
887 }
888 }
889
890 c = cperm[p_col];
891
892 if(r > c) /* Forest Tomlin update */
893 {
894 /* update permutations
895 */
896 j = rorig[c];
897
898 // memmove is more efficient than a for loop
899 // for ( i = c; i < r; ++i )
900 // rorig[i] = rorig[i + 1];
901 memmove(&rorig[c], &rorig[c + 1], (unsigned int)(r - c) * sizeof(int));
902
903 rorig[r] = j;
904
905 for(i = c; i <= r; ++i)
906 rperm[rorig[i]] = i;
907
908 j = corig[c];
909
910 // memmove is more efficient than a for loop
911 // for ( i = c; i < r; ++i )
912 // corig[i] = corig[i + 1];
913 memmove(&corig[c], &corig[c + 1], (unsigned int)(r - c) * sizeof(int));
914
915 corig[r] = j;
916
917 for(i = c; i <= r; ++i)
918 cperm[corig[i]] = i;
919
920
921 rowno = rorig[r];
922
923 j = rbeg[rowno];
924
925 i = rlen[rowno];
926
927 nzCnt -= i;
928
929 if(i < verySparseFactorRat * (dim - c)) // few nonzeros to be eliminated
930 {
931 /**
932 * The following assert is obviously violated if this method is called
933 * with nonzero==0.
934 *
935 * @todo Use an extra member variable as a buffer for the heap instead of
936 * misusing nonz and num. The method enQueueMinRat() seems to
937 * sort the nonzeros or something, for which it only needs
938 * some empty vector of size num.
939 */
940 assert(nonz != 0);
941
942 /* move row r from U to p_work
943 */
944 num = 0;
945
946 for(i += j - 1; i >= j; --i)
947 {
948 k = ridx[i];
949 p_work[k] = rval[i];
950 enQueueMinRat(nonz, &num, cperm[k]);
951 m = --(clen[k]) + cbeg[k];
952
953 for(h = m; cidx[h] != rowno; --h)
954 ;
955
956 cidx[h] = cidx[m];
957
958 cval[h] = cval[m];
959 }
960
961
962 /* Eliminate row r from U to L file
963 */
964 ll = makeLvec(r - c, rowno);
965
966 VectorRational& lval = l.val;
967 int* lidx = l.idx;
968
969 assert((num == 0) || (nonz != 0));
970
971 /* for(i = c; i < r; ++i) */
972 while(num)
973 {
974 #ifndef NDEBUG
975 // The numbers seem to be often 1e-100, is this ok ?
976
977 for(i = 0; i < num; ++i)
978 assert(p_work[corig[nonz[i]]] != 0);
979
980 #endif // NDEBUG
981 i = deQueueMinRat(nonz, &num);
982
983 if(i == r)
984 break;
985
986 k = corig[i];
987
988 assert(p_work[k] != 0);
989
990 n = rorig[i];
991
992 x = p_work[k] * diag[n];
993
994 lidx[ll] = n;
995
996 lval[ll] = x;
997
998 p_work[k] = 0;
999
1000 ll++;
1001
1002 if(spxAbs(x) > l_maxabs)
1003 l_maxabs = spxAbs(x);
1004
1005 j = rbeg[n];
1006
1007 m = rlen[n] + j;
1008
1009 for(; j < m; ++j)
1010 {
1011 int jj = ridx[j];
1012 Rational y = p_work[jj];
1013
1014 if(y == 0)
1015 enQueueMinRat(nonz, &num, cperm[jj]);
1016
1017 y -= x * rval[j];
1018
1019 //p_work[jj] = y + (( y == 0 ) ? SOPLEX_MARKER : 0 );
1020 p_work[jj] = y;
1021 }
1022 }
1023
1024 if(lbeg[l.firstUnused - 1] == ll)
1025 (l.firstUnused)--;
1026 else
1027 lbeg[l.firstUnused] = ll;
1028
1029
1030 /* Set diagonal value
1031 */
1032 if(i != r)
1033 {
1034 stat = SLinSolverRational::SINGULAR;
1035 throw SPxStatusException("XFORE01 The loaded matrix is singular");
1036 }
1037
1038 k = corig[r];
1039
1040 x = p_work[k];
1041 diag[rowno] = 1 / x;
1042 p_work[k] = 0;
1043
1044
1045 /* make row large enough to fit all nonzeros.
1046 */
1047
1048 if(rmax[rowno] < num)
1049 {
1050 rlen[rowno] = 0;
1051 remaxRow(rowno, num);
1052 rval = u.row.val;
1053 ridx = u.row.idx;
1054 }
1055
1056 nzCnt += num;
1057
1058 /* Insert work to updated row thereby clearing work;
1059 */
1060 n = rbeg[rowno];
1061
1062 for(i = 0; i < num; ++i)
1063 {
1064 j = corig[nonz[i]];
1065 x = p_work[j];
1066
1067 // BH 2005-08-24: This if is very important. It may well happen that
1068 // during the elimination of row r a nonzero elements cancels out
1069 // and becomes zero. This would lead to an infinite loop in the
1070 // above elimination code, since the corresponding column index would
1071 // be enqueued for further elimination again and agian.
1072
1073 if(x != 0)
1074 {
1075 if(spxAbs(x) > l_maxabs)
1076 l_maxabs = spxAbs(x);
1077
1078 ridx[n] = j;
1079
1080 rval[n] = x;
1081
1082 p_work[j] = 0;
1083
1084 ++n;
1085
1086 if(clen[j] >= cmax[j])
1087 {
1088 forestReMaxCol(j, clen[j] + 1);
1089 cidx = u.col.idx;
1090 cval = u.col.val;
1091 }
1092
1093 cval[cbeg[j] + clen[j]] = x;
1094
1095 cidx[cbeg[j] + clen[j]++] = rowno;
1096 }
1097 }
1098
1099 rlen[rowno] = n - rbeg[rowno];
1100 }
1101 else /* few nonzeros to be eliminated */
1102 {
1103 /* move row r from U to p_work
1104 */
1105 for(i += j - 1; i >= j; --i)
1106 {
1107 k = ridx[i];
1108 p_work[k] = rval[i];
1109 m = --(clen[k]) + cbeg[k];
1110
1111 for(h = m; cidx[h] != rowno; --h)
1112 ;
1113
1114 cidx[h] = cidx[m];
1115
1116 cval[h] = cval[m];
1117 }
1118
1119
1120 /* Eliminate row r from U to L file
1121 */
1122 ll = makeLvec(r - c, rowno);
1123
1124 VectorRational& lval = l.val;
1125 int* lidx = l.idx;
1126
1127 for(i = c; i < r; ++i)
1128 {
1129 k = corig[i];
1130
1131 if(p_work[k] != 0)
1132 {
1133 n = rorig[i];
1134 x = p_work[k] * diag[n];
1135 lidx[ll] = n;
1136 lval[ll] = x;
1137 p_work[k] = 0;
1138 ll++;
1139
1140 if(spxAbs(x) > l_maxabs)
1141 l_maxabs = spxAbs(x);
1142
1143 j = rbeg[n];
1144
1145 m = rlen[n] + j;
1146
1147 for(; j < m; ++j)
1148 p_work[ridx[j]] -= x * rval[j];
1149 }
1150 }
1151
1152 if(lbeg[l.firstUnused - 1] == ll)
1153 (l.firstUnused)--;
1154 else
1155 lbeg[l.firstUnused] = ll;
1156
1157
1158 /* Set diagonal value
1159 */
1160 k = corig[r];
1161
1162 x = p_work[k];
1163
1164 if(x == 0)
1165 {
1166 stat = SLinSolverRational::SINGULAR;
1167 throw SPxStatusException("XFORE02 The loaded matrix is singular");
1168 // return;
1169 }
1170
1171 diag[rowno] = 1 / x;
1172
1173 p_work[k] = 0;
1174
1175
1176 /* count remaining nonzeros in work and make row large enough
1177 * to fit them all.
1178 */
1179 n = 0;
1180
1181 for(i = r + 1; i < dim; ++i)
1182 if(p_work[corig[i]] != 0)
1183 n++;
1184
1185 if(rmax[rowno] < n)
1186 {
1187 rlen[rowno] = 0;
1188 remaxRow(rowno, n);
1189 rval = u.row.val;
1190 ridx = u.row.idx;
1191 }
1192
1193 nzCnt += n;
1194
1195 /* Insert p_work to updated row thereby clearing p_work;
1196 */
1197 n = rbeg[rowno];
1198
1199 for(i = r + 1; i < dim; ++i)
1200 {
1201 j = corig[i];
1202 x = p_work[j];
1203
1204 if(x != 0)
1205 {
1206 if(spxAbs(x) > l_maxabs)
1207 l_maxabs = spxAbs(x);
1208
1209 ridx[n] = j;
1210
1211 rval[n] = x;
1212
1213 p_work[j] = 0;
1214
1215 ++n;
1216
1217 if(clen[j] >= cmax[j])
1218 {
1219 forestReMaxCol(j, clen[j] + 1);
1220 cidx = u.col.idx;
1221 cval = u.col.val;
1222 }
1223
1224 cval[cbeg[j] + clen[j]] = x;
1225
1226 cidx[cbeg[j] + clen[j]++] = rowno;
1227 }
1228 }
1229
1230 rlen[rowno] = n - rbeg[rowno];
1231 }
1232 }
1233
1234 else if(r == c)
1235 {
1236 /* Move diagonal element to diag. Note, that it must be the last
1237 * element, since it has just been inserted above.
1238 */
1239 rowno = rorig[r];
1240 i = rbeg[rowno] + --(rlen[rowno]);
1241 diag[rowno] = 1 / rval[i];
1242
1243 for(j = i = --(clen[p_col]) + cbeg[p_col]; cidx[i] != rowno; --i)
1244 ;
1245
1246 cidx[i] = cidx[j];
1247
1248 cval[i] = cval[j];
1249 }
1250 else /* r < c */
1251 {
1252 stat = SLinSolverRational::SINGULAR;
1253 throw SPxStatusException("XFORE03 The loaded matrix is singular");
1254 // return;
1255 }
1256
1257 maxabs = l_maxabs;
1258
1259 assert(isConsistent());
1260 stat = SLinSolverRational::OK;
1261 }
1262
1263 inline void CLUFactorRational::update(int p_col, Rational* p_work, const int* p_idx, int num)
1264 {
1265 int ll, i, j;
1266 Rational x, rezi;
1267
1268 assert(p_work[p_col] != 0);
1269 rezi = 1 / p_work[p_col];
1270 p_work[p_col] = 0;
1271
1272 ll = makeLvec(num, p_col);
1273 // ll = fac->makeLvec(num, col);
1274 VectorRational& lval = l.val;
1275 int* lidx = l.idx;
1276
1277 for(i = num - 1; (j = p_idx[i]) != p_col; --i)
1278 {
1279 lidx[ll] = j;
1280 lval[ll] = rezi * p_work[j];
1281 p_work[j] = 0;
1282 ++ll;
1283 }
1284
1285 lidx[ll] = p_col;
1286
1287 lval[ll] = 1 - rezi;
1288 ++ll;
1289
1290 for(--i; i >= 0; --i)
1291 {
1292 j = p_idx[i];
1293 lidx[ll] = j;
1294 lval[ll] = x = rezi * p_work[j];
1295 p_work[j] = 0;
1296 ++ll;
1297
1298 if(spxAbs(x) > maxabs)
1299 maxabs = spxAbs(x);
1300 }
1301
1302 stat = SLinSolverRational::OK;
1303 }
1304
1305 inline void CLUFactorRational::updateNoClear(
1306 int p_col,
1307 const Rational* p_work,
1308 const int* p_idx,
1309 int num)
1310 {
1311 int ll, i, j;
1312 Rational x, rezi;
1313
1314 assert(p_work[p_col] != 0);
1315 rezi = 1 / p_work[p_col];
1316 ll = makeLvec(num, p_col);
1317 //ll = fac->makeLvec(num, col);
1318 VectorRational& lval = l.val;
1319 int* lidx = l.idx;
1320
1321 for(i = num - 1; (j = p_idx[i]) != p_col; --i)
1322 {
1323 lidx[ll] = j;
1324 lval[ll] = rezi * p_work[j];
1325 ++ll;
1326 }
1327
1328 lidx[ll] = p_col;
1329
1330 lval[ll] = 1 - rezi;
1331 ++ll;
1332
1333 for(--i; i >= 0; --i)
1334 {
1335 j = p_idx[i];
1336 lidx[ll] = j;
1337 lval[ll] = x = rezi * p_work[j];
1338 ++ll;
1339
1340 if(spxAbs(x) > maxabs)
1341 maxabs = spxAbs(x);
1342 }
1343
1344 stat = SLinSolverRational::OK;
1345 }
1346
1347 /*****************************************************************************/
1348 /*
1349 * Temporary data structures.
1350 */
1351
1352 /*
1353 * For the i=th column the situation might look like this:
1354 *
1355 * \begin{verbatim}
1356 * idx = ....................iiiIIIIII-----..............
1357 * cbeg[i] = ^
1358 * cact[i] = +----+
1359 * clen[i] = +-------+
1360 * cmax[i] = +------------+
1361 *
1362 * Indices clen[i]-cbeg[i]: ^^^
1363 * \end{verbatim}
1364 * belong to column i, but have already been pivotal and don't belong to
1365 * the active submatrix.
1366 */
1367
1368 /****************************************************************************/
1369 /*
1370 * Initialize row and column file of working matrix and
1371 * mark column singletons.
1372 */
1373 inline void CLUFactorRational::initFactorMatrix(const SVectorRational** vec)
1374 {
1375
1376 Rational x;
1377 int m;
1378 int tot;
1379 Dring* rring, *lastrring;
1380 Dring* cring, *lastcring;
1381 const SVectorRational* psv;
1382 int* sing = temp.s_mark;
1383
1384 /* Initialize:
1385 * - column file thereby remembering column singletons in |sing|.
1386 * - nonzeros counts per row
1387 * - total number of nonzeros
1388 */
1389
1390 for(int i = 0; i < thedim; i++)
1391 u.row.max[i] = u.row.len[i] = 0;
1392
1393 tot = 0;
1394
1395 for(int i = 0; i < thedim; i++)
1396 {
1397 int k;
1398
1399 psv = vec[i];
1400 k = psv->size();
1401
1402 if(k > 1)
1403 {
1404 tot += k;
1405
1406 for(int j = 0; j < k; ++j)
1407 u.row.max[psv->index(j)]++;
1408 }
1409 else if(k == 0)
1410 {
1411 stat = SLinSolverRational::SINGULAR;
1412 return;
1413 }
1414 }
1415
1416 /* Resize nonzero memory if necessary
1417 */
1418 minRowMem(int(rowMemMult * tot));
1419
1420 minColMem(int(colMemMult * tot));
1421
1422 minLMem(int(lMemMult * tot));
1423
1424
1425 /* Initialize:
1426 * - row ring lists
1427 * - row vectors in file
1428 * - column ring lists
1429 */
1430 u.row.start[0] = 0;
1431
1432 rring = u.row.elem;
1433
1434 lastrring = &(u.row.list);
1435
1436 lastrring->idx = thedim;
1437
1438 lastrring->next = rring;
1439
1440 cring = u.col.elem;
1441
1442 lastcring = &(u.col.list);
1443
1444 lastcring->idx = thedim;
1445
1446 lastcring->next = cring;
1447
1448 m = 0;
1449
1450 for(int i = 0; i < thedim; i++)
1451 {
1452 u.row.start[i] = m;
1453 m += u.row.max[i];
1454
1455 rring->idx = i;
1456 rring->prev = lastrring;
1457 lastrring->next = rring;
1458 lastrring = rring;
1459 ++rring;
1460
1461 cring->idx = i;
1462 cring->prev = lastcring;
1463 lastcring->next = cring;
1464 lastcring = cring;
1465 ++cring;
1466 }
1467
1468 u.row.start[thedim] = 0;
1469
1470 u.row.max[thedim] = 0;
1471 u.row.used = m;
1472
1473 lastrring->next = &(u.row.list);
1474 lastrring->next->prev = lastrring;
1475
1476 lastcring->next = &(u.col.list);
1477 lastcring->next->prev = lastcring;
1478
1479 /* Copy matrix to row and column file
1480 * excluding and marking column singletons!
1481 */
1482 m = 0;
1483 temp.stage = 0;
1484
1485 initMaxabs = 0;
1486
1487 for(int i = 0; i < thedim; i++)
1488 {
1489 int nnonzeros;
1490
1491 psv = vec[i];
1492 u.col.start[i] = m;
1493
1494 /* check whether number of nonzeros above tolerance is 0, 1 or >= 2 */
1495 nnonzeros = 0;
1496
1497 for(int j = 0; j < psv->size() && nnonzeros <= 1; j++)
1498 {
1499 if(psv->value(j) != 0)
1500 nnonzeros++;
1501 }
1502
1503 /* basis is singular due to empty column */
1504 if(nnonzeros == 0)
1505 {
1506 stat = SLinSolverRational::SINGULAR;
1507 return;
1508 }
1509
1510 /* exclude column singletons */
1511 else if(nnonzeros == 1)
1512 {
1513 int j;
1514
1515 /* find nonzero */
1516
1517 for(j = 0; psv->value(j) == 0; j++)
1518 ;
1519
1520 assert(j < psv->size());
1521
1522 /* basis is singular due to two linearly dependent column singletons */
1523 if(row.perm[psv->index(j)] >= 0)
1524 {
1525 stat = SLinSolverRational::SINGULAR;
1526 return;
1527 }
1528
1529 /* update maximum absolute nonzero value */
1530 x = psv->value(j);
1531
1532 if(spxAbs(x) > initMaxabs)
1533 initMaxabs = spxAbs(x);
1534
1535 /* permute to front and mark as singleton */
1536 setPivot(temp.stage, i, psv->index(j), x);
1537
1538 sing[temp.stage] = i;
1539
1540 temp.stage++;
1541
1542 /* set column length to zero */
1543 temp.s_cact[i] = u.col.len[i] = u.col.max[i] = 0;
1544 }
1545
1546 /* add to active matrix if not a column singleton */
1547 else
1548 {
1549 int end;
1550 int k;
1551
1552 /* go through all nonzeros in column */
1553 assert(nnonzeros >= 2);
1554 nnonzeros = 0;
1555
1556 for(int j = 0; j < psv->size(); j++)
1557 {
1558 x = psv->value(j);
1559
1560 if(x != 0)
1561 {
1562 /* add to column array */
1563 k = psv->index(j);
1564 u.col.idx[m] = k;
1565 m++;
1566
1567 /* add to row array */
1568 end = u.row.start[k] + u.row.len[k];
1569 u.row.idx[end] = i;
1570 u.row.val[end] = x;
1571 u.row.len[k]++;
1572
1573 /* update maximum absolute nonzero value */
1574
1575 if(spxAbs(x) > initMaxabs)
1576 initMaxabs = spxAbs(x);
1577
1578 nnonzeros++;
1579 }
1580 }
1581
1582 assert(nnonzeros >= 2);
1583
1584 /* set column length */
1585 temp.s_cact[i] = u.col.len[i] = u.col.max[i] = nnonzeros;
1586 }
1587 }
1588
1589 u.col.used = m;
1590 }
1591
1592
1593
1594 /*****************************************************************************/
1595 /*
1596 * Remove column singletons
1597 */
1598
1599 inline void CLUFactorRational::colSingletons()
1600 {
1601 int i, j, k, n;
1602 int len;
1603 int p_col, p_row, newrow;
1604 int* idx;
1605 int* rorig = row.orig;
1606 int* rperm = row.perm;
1607 int* sing = temp.s_mark;
1608
1609
1610 /* Iteratively update column counts due to removed column singletons
1611 * thereby removing new arising columns singletons
1612 * and computing the index of the first row singleton (-1)
1613 * until no more can be found.
1614 */
1615
1616 for(i = 0; i < temp.stage; ++i)
1617 {
1618 p_row = rorig[i];
1619 assert(p_row >= 0);
1620 idx = &(u.row.idx[u.row.start[p_row]]);
1621 len = u.row.len[p_row];
1622
1623 for(j = 0; j < len; ++j)
1624 {
1625 /* Move pivotal nonzeros to front of column.
1626 */
1627 p_col = idx[j];
1628 assert(temp.s_cact[p_col] > 0);
1629
1630 n = u.col.start[p_col] + u.col.len[p_col] - temp.s_cact[p_col];
1631
1632 for(k = n; u.col.idx[k] != p_row; ++k)
1633 ;
1634
1635 assert(k < u.col.start[p_col] + u.col.len[p_col]);
1636
1637 u.col.idx[k] = u.col.idx[n];
1638
1639 u.col.idx[n] = p_row;
1640
1641 n = --(temp.s_cact[p_col]); /* column nonzeros of ACTIVE matrix */
1642
1643 if(n == 1) /* Here is another singleton */
1644 {
1645 newrow = u.col.idx[--u.col.len[p_col] + u.col.start[p_col]];
1646
1647 /* Ensure, matrix not singular
1648 */
1649
1650 if(rperm[newrow] >= 0)
1651 {
1652 stat = SLinSolverRational::SINGULAR;
1653 return;
1654 }
1655
1656 /* Find singleton in row.
1657 */
1658 n = u.row.start[newrow] + (--(u.row.len[newrow]));
1659
1660 for(k = n; u.row.idx[k] != p_col; --k)
1661 ;
1662
1663 /* Remove singleton from column.
1664 */
1665 setPivot(temp.stage, p_col, newrow, u.row.val[k]);
1666
1667 sing[temp.stage++] = p_col;
1668
1669 /* Move pivot element to diag.
1670 */
1671 u.row.val[k] = u.row.val[n];
1672
1673 u.row.idx[k] = u.row.idx[n];
1674 }
1675 else if(n == 0)
1676 {
1677 stat = SLinSolverRational::SINGULAR;
1678 return;
1679 }
1680 }
1681 }
1682
1683 assert(temp.stage <= thedim);
1684 }
1685
1686
1687 /*****************************************************************************/
1688 /*
1689 * Remove row singletons
1690 */
1691 inline void CLUFactorRational::rowSingletons()
1692 {
1693 Rational pval;
1694 int i, j, k, ll, r;
1695 int p_row, p_col, len, rs, lk;
1696 int* idx;
1697 int* rperm = row.perm;
1698 int* sing = temp.s_mark;
1699
1700 /* Mark row singletons
1701 */
1702 rs = temp.stage;
1703
1704 for(i = 0; i < thedim; ++i)
1705 {
1706 if(rperm[i] < 0 && u.row.len[i] == 1)
1707 sing[temp.stage++] = i;
1708 }
1709
1710 /* Eliminate row singletons
1711 * thereby marking newly arising ones
1712 * until no more can be found.
1713 */
1714 for(; rs < temp.stage; ++rs)
1715 {
1716 /* Move pivot element from row file to diag
1717 */
1718 p_row = sing[rs];
1719 j = u.row.start[p_row];
1720 p_col = u.row.idx[j];
1721 pval = u.row.val[j];
1722 setPivot(rs, p_col, p_row, pval);
1723 u.row.len[p_row] = 0;
1724
1725 /* Remove pivot column form workingmatrix
1726 * thereby building up L vector.
1727 */
1728 idx = &(u.col.idx[u.col.start[p_col]]);
1729 i = temp.s_cact[p_col]; /* nr. nonzeros of new L vector */
1730 lk = makeLvec(i - 1, p_row);
1731 len = u.col.len[p_col];
1732 i = (u.col.len[p_col] -= i); /* remove pivot column from U */
1733
1734 for(; i < len; ++i)
1735 {
1736 r = idx[i];
1737
1738 if(r != p_row)
1739 {
1740 /* Find pivot column in row.
1741 */
1742 ll = --(u.row.len[r]);
1743 k = u.row.start[r] + ll;
1744
1745 for(j = k; u.row.idx[j] != p_col; --j)
1746 ;
1747
1748 assert(k >= u.row.start[r]);
1749
1750 /* Initialize L vector
1751 */
1752 l.idx[lk] = r;
1753
1754 l.val[lk] = u.row.val[j] / pval;
1755
1756 ++lk;
1757
1758 /* Remove pivot column from row.
1759 */
1760 u.row.idx[j] = u.row.idx[k];
1761
1762 u.row.val[j] = u.row.val[k];
1763
1764 /* Check new row length.
1765 */
1766 if(ll == 1)
1767 sing[temp.stage++] = r;
1768 else if(ll == 0)
1769 {
1770 stat = SLinSolverRational::SINGULAR;
1771 return;
1772 }
1773 }
1774 }
1775 }
1776 }
1777
1778
1779 /*****************************************************************************/
1780 /*
1781 * Init nonzero number Ring lists
1782 * and required entries of arrays max and mark
1783 */
1784
1785 inline void CLUFactorRational::initFactorRings()
1786 {
1787 int i;
1788 int* rperm = row.perm;
1789 int* cperm = col.perm;
1790 CLUFactorRational::Pring* ring;
1791
1792 assert(thedim >= 0);
1793 spx_alloc(temp.pivot_col, thedim + 1);
1794 spx_alloc(temp.pivot_colNZ, thedim + 1);
1795 spx_alloc(temp.pivot_row, thedim + 1);
1796 spx_alloc(temp.pivot_rowNZ, thedim + 1);
1797
1798 for(i = thedim - temp.stage; i >= 0; --i)
1799 {
1800 initDR(temp.pivot_colNZ[i]);
1801 initDR(temp.pivot_rowNZ[i]);
1802 }
1803
1804 for(i = 0; i < thedim; ++i)
1805 {
1806 if(rperm[i] < 0)
1807 {
1808 if(u.row.len[i] <= 0)
1809 {
1810 stat = SLinSolverRational::SINGULAR;
1811 return;
1812 }
1813
1814 ring = &(temp.pivot_rowNZ[u.row.len[i]]);
1815
1816 init2DR(temp.pivot_row[i], *ring);
1817 temp.pivot_row[i].idx = i;
1818 temp.s_max[i] = -1;
1819 }
1820
1821 if(cperm[i] < 0)
1822 {
1823 if(temp.s_cact[i] <= 0)
1824 {
1825 stat = SLinSolverRational::SINGULAR;
1826 return;
1827 }
1828
1829 ring = &(temp.pivot_colNZ[temp.s_cact[i]]);
1830
1831 init2DR(temp.pivot_col[i], *ring);
1832 temp.pivot_col[i].idx = i;
1833 temp.s_mark[i] = 0;
1834 }
1835 }
1836 }
1837
1838 inline void CLUFactorRational::freeFactorRings(void)
1839 {
1840
1841 if(temp.pivot_col)
1842 spx_free(temp.pivot_col);
1843
1844 if(temp.pivot_colNZ)
1845 spx_free(temp.pivot_colNZ);
1846
1847 if(temp.pivot_row)
1848 spx_free(temp.pivot_row);
1849
1850 if(temp.pivot_rowNZ)
1851 spx_free(temp.pivot_rowNZ);
1852 }
1853
1854
1855 /*
1856 * Eliminate all row singletons from nucleus.
1857 * A row singleton may well be column singleton at the same time!
1858 */
1859 inline void CLUFactorRational::eliminateRowSingletons()
1860 {
1861 int i, j, k, ll, r;
1862 int len, lk;
1863 int pcol, prow;
1864 Rational pval;
1865 int* idx;
1866 CLUFactorRational::Pring* sing;
1867
1868 for(sing = temp.pivot_rowNZ[1].prev; sing != &(temp.pivot_rowNZ[1]); sing = sing->prev)
1869 {
1870 prow = sing->idx;
1871 i = u.row.start[prow];
1872 pcol = u.row.idx[i];
1873 pval = u.row.val[i];
1874 setPivot(temp.stage++, pcol, prow, pval);
1875 u.row.len[prow] = 0;
1876 removeDR(temp.pivot_col[pcol]);
1877
1878 /* Eliminate pivot column and build L vector.
1879 */
1880 i = temp.s_cact[pcol];
1881
1882 if(i > 1)
1883 {
1884 idx = &(u.col.idx[u.col.start[pcol]]);
1885 len = u.col.len[pcol];
1886 lk = makeLvec(i - 1, prow);
1887 i = u.col.len[pcol] -= i;
1888
1889 for(; (r = idx[i]) != prow; ++i)
1890 {
1891 /* Find pivot column in row.
1892 */
1893 ll = --(u.row.len[r]);
1894 k = u.row.start[r] + ll;
1895
1896 for(j = k; u.row.idx[j] != pcol; --j)
1897 ;
1898
1899 assert(j >= u.row.start[r]);
1900
1901 /* Initialize L vector
1902 */
1903 l.idx[lk] = r;
1904
1905 l.val[lk] = u.row.val[j] / pval;
1906
1907 ++lk;
1908
1909 /* Remove pivot column from row.
1910 */
1911 u.row.idx[j] = u.row.idx[k];
1912
1913 u.row.val[j] = u.row.val[k];
1914
1915 /* Move column to appropriate nonzero ring.
1916 */
1917 removeDR(temp.pivot_row[r]);
1918
1919 init2DR(temp.pivot_row[r], temp.pivot_rowNZ[ll]);
1920
1921 assert(row.perm[r] < 0);
1922
1923 temp.s_max[r] = -1;
1924 }
1925
1926 /* skip pivot element */
1927 assert(i < len && "ERROR: pivot column does not contain pivot row");
1928
1929 for(++i; i < len; ++i)
1930 {
1931 /* Find pivot column in row.
1932 */
1933 r = idx[i];
1934 ll = --(u.row.len[r]);
1935 k = u.row.start[r] + ll;
1936
1937 for(j = k; u.row.idx[j] != pcol; --j)
1938 ;
1939
1940 assert(j >= u.row.start[r]);
1941
1942 /* Initialize L vector
1943 */
1944 l.idx[lk] = r;
1945
1946 l.val[lk] = u.row.val[j] / pval;
1947
1948 ++lk;
1949
1950 /* Remove pivot column from row.
1951 */
1952 u.row.idx[j] = u.row.idx[k];
1953
1954 u.row.val[j] = u.row.val[k];
1955
1956 /* Move column to appropriate nonzero ring.
1957 */
1958 removeDR(temp.pivot_row[r]);
1959
1960 init2DR(temp.pivot_row[r], temp.pivot_rowNZ[ll]);
1961
1962 assert(row.perm[r] < 0);
1963
1964 temp.s_max[r] = -1;
1965 }
1966 }
1967 else
1968 u.col.len[pcol] -= i;
1969 }
1970
1971 initDR(temp.pivot_rowNZ[1]); /* Remove all row singletons from list */
1972 }
1973
1974
1975
1976 /*
1977 * Eliminate all column singletons from nucleus.
1978 * A column singleton must not be row singleton at the same time!
1979 */
1980 inline void CLUFactorRational::eliminateColSingletons()
1981 {
1982 int i, j, k, m, c;
1983 int pcol, prow;
1984 CLUFactorRational::Pring* sing;
1985
1986 for(sing = temp.pivot_colNZ[1].prev;
1987 sing != &(temp.pivot_colNZ[1]);
1988 sing = sing->prev)
1989 {
1990 /* Find pivot value
1991 */
1992 pcol = sing->idx;
1993 j = --(u.col.len[pcol]) + u.col.start[pcol]; /* remove pivot column */
1994 prow = u.col.idx[j];
1995 removeDR(temp.pivot_row[prow]);
1996
1997 j = --(u.row.len[prow]) + u.row.start[prow];
1998
1999 for(i = j; (c = u.row.idx[i]) != pcol; --i)
2000 {
2001 m = u.col.len[c] + u.col.start[c] - (temp.s_cact[c])--;
2002
2003 for(k = m; u.col.idx[k] != prow; ++k)
2004 ;
2005
2006 u.col.idx[k] = u.col.idx[m];
2007
2008 u.col.idx[m] = prow;
2009
2010 m = temp.s_cact[c];
2011
2012 removeDR(temp.pivot_col[c]);
2013
2014 init2DR(temp.pivot_col[c], temp.pivot_colNZ[m]);
2015
2016 assert(col.perm[c] < 0);
2017 }
2018
2019 /* remove pivot element from pivot row
2020 */
2021 setPivot(temp.stage++, pcol, prow, u.row.val[i]);
2022
2023 u.row.idx[i] = u.row.idx[j];
2024
2025 u.row.val[i] = u.row.val[j];
2026
2027 j = u.row.start[prow];
2028
2029 for(--i; i >= j; --i)
2030 {
2031 c = u.row.idx[i];
2032 m = u.col.len[c] + u.col.start[c] - (temp.s_cact[c])--;
2033
2034 for(k = m; u.col.idx[k] != prow; ++k)
2035 ;
2036
2037 u.col.idx[k] = u.col.idx[m];
2038
2039 u.col.idx[m] = prow;
2040
2041 m = temp.s_cact[c];
2042
2043 removeDR(temp.pivot_col[c]);
2044
2045 init2DR(temp.pivot_col[c], temp.pivot_colNZ[m]);
2046
2047 assert(col.perm[c] < 0);
2048 }
2049 }
2050
2051 initDR(temp.pivot_colNZ[1]); /* Remove all column singletons from list */
2052 }
2053
2054 /*
2055 * No singletons available: Select pivot elements.
2056 */
2057 inline void CLUFactorRational::selectPivots(const Rational& threshold)
2058 {
2059 int ii;
2060 int i;
2061 int j;
2062 int k;
2063 int ll = -1; // This value should never be used.
2064 int kk;
2065 int m;
2066 int count;
2067 int num;
2068 int rw = -1; // This value should never be used.
2069 int cl = -1; // This value should never be used.
2070 int len;
2071 int beg;
2072 Rational l_maxabs;
2073 Rational x = 0; // This value should never be used.
2074 int mkwtz;
2075 int candidates;
2076
2077 candidates = thedim - temp.stage - 1;
2078
2079 if(candidates > 4)
2080 candidates = 4;
2081
2082 num = 0;
2083
2084 count = 2;
2085
2086 for(;;)
2087 {
2088 ii = -1;
2089
2090 if(temp.pivot_rowNZ[count].next != &(temp.pivot_rowNZ[count]))
2091 {
2092 rw = temp.pivot_rowNZ[count].next->idx;
2093 beg = u.row.start[rw];
2094 len = u.row.len[rw] + beg - 1;
2095
2096 /* set l_maxabs to maximum absolute value in row
2097 * (compute it if necessary).
2098 */
2099
2100 if(temp.s_max[rw] < 0)
2101 {
2102 l_maxabs = spxAbs(u.row.val[len]);
2103
2104 for(i = len - 1; i >= beg; --i)
2105 if(l_maxabs < spxAbs(u.row.val[i]))
2106 l_maxabs = spxAbs(u.row.val[i]);
2107
2108 temp.s_max[rw] = l_maxabs; /* ##### */
2109 }
2110 else
2111 l_maxabs = temp.s_max[rw];
2112
2113 l_maxabs *= threshold;
2114
2115 /* select pivot element with lowest markowitz number in row
2116 */
2117 mkwtz = thedim + 1;
2118
2119 for(i = len; i >= beg; --i)
2120 {
2121 k = u.row.idx[i];
2122 j = temp.s_cact[k];
2123 x = u.row.val[i];
2124
2125 if(j < mkwtz && spxAbs(x) > l_maxabs)
2126 {
2127 mkwtz = j;
2128 cl = k;
2129 ii = i;
2130
2131 if(j <= count) /* ##### */
2132 break;
2133 }
2134 }
2135 }
2136 else if(temp.pivot_colNZ[count].next != &(temp.pivot_colNZ[count]))
2137 {
2138 cl = temp.pivot_colNZ[count].next->idx;
2139 beg = u.col.start[cl];
2140 len = u.col.len[cl] + beg - 1;
2141 beg = len - temp.s_cact[cl] + 1;
2142 assert(count == temp.s_cact[cl]);
2143
2144 /* select pivot element with lowest markowitz number in column
2145 */
2146 mkwtz = thedim + 1;
2147
2148 for(i = len; i >= beg; --i)
2149 {
2150 k = u.col.idx[i];
2151 j = u.row.len[k];
2152
2153 if(j < mkwtz)
2154 {
2155 /* ensure that element (cl,k) is stable.
2156 */
2157 if(temp.s_max[k] > 0)
2158 {
2159 /* case 1: l_maxabs is known
2160 */
2161 for(m = u.row.start[k], kk = m + u.row.len[k] - 1;
2162 kk >= m; --kk)
2163 {
2164 if(u.row.idx[kk] == cl)
2165 {
2166 x = u.row.val[kk];
2167 ll = kk;
2168 break;
2169 }
2170 }
2171
2172 l_maxabs = temp.s_max[k];
2173 }
2174 else
2175 {
2176 /* case 2: l_maxabs needs to be computed
2177 */
2178 m = u.row.start[k];
2179 l_maxabs = spxAbs(u.row.val[m]);
2180
2181 for(kk = m + u.row.len[k] - 1; kk >= m; --kk)
2182 {
2183 if(l_maxabs < spxAbs(u.row.val[kk]))
2184 l_maxabs = spxAbs(u.row.val[kk]);
2185
2186 if(u.row.idx[kk] == cl)
2187 {
2188 x = u.row.val[kk];
2189 ll = kk;
2190 break;
2191 }
2192 }
2193
2194 for(--kk; kk > m; --kk)
2195 {
2196 if(l_maxabs < spxAbs(u.row.val[kk]))
2197 l_maxabs = spxAbs(u.row.val[kk]);
2198 }
2199
2200 temp.s_max[k] = l_maxabs;
2201 }
2202
2203 l_maxabs *= threshold;
2204
2205 if(spxAbs(x) > l_maxabs)
2206 {
2207 mkwtz = j;
2208 rw = k;
2209 ii = ll;
2210
2211 if(j <= count + 1)
2212 break;
2213 }
2214 }
2215 }
2216 }
2217 else
2218 {
2219 ++count;
2220 continue;
2221 }
2222
2223 assert(cl >= 0);
2224
2225 removeDR(temp.pivot_col[cl]);
2226 initDR(temp.pivot_col[cl]);
2227
2228 if(ii >= 0)
2229 {
2230 /* Initialize selected pivot element
2231 */
2232 CLUFactorRational::Pring* pr;
2233 temp.pivot_row[rw].pos = ii - u.row.start[rw];
2234 temp.pivot_row[rw].mkwtz = mkwtz = (mkwtz - 1) * (count - 1);
2235 // ??? mkwtz originally was long,
2236 // maybe to avoid an overflow in this instruction?
2237
2238 for(pr = temp.pivots.next; pr->idx >= 0; pr = pr->next)
2239 {
2240 if(pr->idx == rw || pr->mkwtz >= mkwtz)
2241 break;
2242 }
2243
2244 pr = pr->prev;
2245
2246 if(pr->idx != rw)
2247 {
2248 removeDR(temp.pivot_row[rw]);
2249 init2DR(temp.pivot_row[rw], *pr);
2250 }
2251
2252 num++;
2253
2254 if(num >= candidates)
2255 break;
2256 }
2257 }
2258
2259 /*
2260 * while(temp.temp.next->mkwtz < temp.temp.prev->mkwtz)
2261 * {
2262 * Pring *pr;
2263 * pr = temp.temp.prev;
2264 * removeDR(*pr);
2265 * init2DR (*pr, rowNZ[u.row.len[pr->idx]]);
2266 }
2267 */
2268
2269 assert(row.perm[rw] < 0);
2270
2271 assert(col.perm[cl] < 0);
2272 }
2273
2274
2275 /*
2276 * Perform L and update loop for row r
2277 */
2278 inline int CLUFactorRational::updateRow(int r,
2279 int lv,
2280 int prow,
2281 int pcol,
2282 const Rational& pval)
2283 {
2284 int fill;
2285 Rational x, lx;
2286 int c, i, j, k, ll, m, n;
2287
2288 n = u.row.start[r];
2289 m = --(u.row.len[r]) + n;
2290
2291 /* compute L vector entry and
2292 * and remove pivot column form row file
2293 */
2294
2295 for(j = m; u.row.idx[j] != pcol; --j)
2296 ;
2297
2298 lx = u.row.val[j] / pval;
2299
2300 l.val[lv] = lx;
2301
2302 l.idx[lv] = r;
2303
2304 ++lv;
2305
2306 u.row.idx[j] = u.row.idx[m];
2307
2308 u.row.val[j] = u.row.val[m];
2309
2310
2311 /* update loop (I) and
2312 * computing expected fill
2313 */
2314 fill = u.row.len[prow];
2315
2316 for(j = m - 1; j >= n; --j)
2317 {
2318 c = u.row.idx[j];
2319
2320 if(temp.s_mark[c])
2321 {
2322 /* count fill elements.
2323 */
2324 temp.s_mark[c] = 0;
2325 --fill;
2326
2327 /* update row values
2328 */
2329 x = u.row.val[j] -= work[c] * lx;
2330
2331 if(x == 0)
2332 {
2333 /* Eliminate zero from row r
2334 */
2335 --u.row.len[r];
2336 --m;
2337 u.row.val[j] = u.row.val[m];
2338 u.row.idx[j] = u.row.idx[m];
2339
2340 /* Eliminate zero from column c
2341 */
2342 --(temp.s_cact[c]);
2343 k = --(u.col.len[c]) + u.col.start[c];
2344
2345 for(i = k; u.col.idx[i] != r; --i)
2346 ;
2347
2348 u.col.idx[i] = u.col.idx[k];
2349 }
2350 }
2351 }
2352
2353
2354 /* create space for fill in row file
2355 */
2356 ll = u.row.len[r];
2357
2358 if(ll + fill > u.row.max[r])
2359 remaxRow(r, ll + fill);
2360
2361 ll += u.row.start[r];
2362
2363 /* fill creating update loop (II)
2364 */
2365 for(j = u.row.start[prow], m = j + u.row.len[prow]; j < m; ++j)
2366 {
2367 c = u.row.idx[j];
2368
2369 if(temp.s_mark[c])
2370 {
2371 x = - work[c] * lx;
2372
2373 if(x != 0)
2374 {
2375 /* produce fill element in row r
2376 */
2377 u.row.val[ll] = x;
2378 u.row.idx[ll] = c;
2379 ll++;
2380 u.row.len[r]++;
2381
2382 /* produce fill element in column c
2383 */
2384
2385 if(u.col.len[c] >= u.col.max[c])
2386 remaxCol(c, u.col.len[c] + 1);
2387
2388 u.col.idx[u.col.start[c] + (u.col.len[c])++] = r;
2389
2390 temp.s_cact[c]++;
2391 }
2392 }
2393 else
2394 temp.s_mark[c] = 1;
2395 }
2396
2397 /* move row to appropriate list.
2398 */
2399 removeDR(temp.pivot_row[r]);
2400
2401 init2DR(temp.pivot_row[r], temp.pivot_rowNZ[u.row.len[r]]);
2402
2403 assert(row.perm[r] < 0);
2404
2405 temp.s_max[r] = -1;
2406
2407 return lv;
2408 }
2409
2410 /*
2411 * Eliminate pivot element
2412 */
2413 inline void CLUFactorRational::eliminatePivot(int prow, int pos)
2414 {
2415 int i, j, k, m = -1;
2416 int lv = -1; // This value should never be used.
2417 int pcol;
2418 Rational pval;
2419 int pbeg = u.row.start[prow];
2420 int plen = --(u.row.len[prow]);
2421 int pend = pbeg + plen;
2422
2423
2424 /* extract pivot element */
2425 i = pbeg + pos;
2426 pcol = u.row.idx[i];
2427 pval = u.row.val[i];
2428 removeDR(temp.pivot_col[pcol]);
2429 initDR(temp.pivot_col[pcol]);
2430
2431 /* remove pivot from pivot row */
2432 u.row.idx[i] = u.row.idx[pend];
2433 u.row.val[i] = u.row.val[pend];
2434
2435 /* set pivot element and construct L vector */
2436 setPivot(temp.stage++, pcol, prow, pval);
2437
2438 /**@todo If this test failes, lv has no value. I suppose that in this
2439 * case none of the loops below that uses lv is executed.
2440 * But this is unproven.
2441 */
2442
2443 if(temp.s_cact[pcol] - 1 > 0)
2444 lv = makeLvec(temp.s_cact[pcol] - 1, prow);
2445
2446 /* init working vector,
2447 * remove pivot row from working matrix
2448 * and remove columns from list.
2449 */
2450 for(i = pbeg; i < pend; ++i)
2451 {
2452 j = u.row.idx[i];
2453 temp.s_mark[j] = 1;
2454 work[j] = u.row.val[i];
2455 removeDR(temp.pivot_col[j]);
2456 m = u.col.start[j] + u.col.len[j] - temp.s_cact[j];
2457
2458 for(k = m; u.col.idx[k] != prow; ++k)
2459 ;
2460
2461 u.col.idx[k] = u.col.idx[m];
2462
2463 u.col.idx[m] = prow;
2464
2465 temp.s_cact[j]--;
2466 }
2467
2468 /* perform L and update loop
2469 */
2470 for(i = u.col.len[pcol] - temp.s_cact[pcol];
2471 (m = u.col.idx[u.col.start[pcol] + i]) != prow;
2472 ++i)
2473 {
2474 assert(row.perm[m] < 0);
2475 assert(lv >= 0);
2476 /* coverity[negative_returns] */
2477 updateRow(m, lv++, prow, pcol, pval);
2478 }
2479
2480 /* skip pivot row */
2481
2482 m = u.col.len[pcol];
2483
2484 for(++i; i < m; ++i)
2485 {
2486 assert(lv >= 0);
2487 /* coverity[negative_returns] */
2488 updateRow(u.col.idx[u.col.start[pcol] + i], lv++, prow, pcol, pval);
2489 }
2490
2491 /* remove pivot column from column file.
2492 */
2493 u.col.len[pcol] -= temp.s_cact[pcol];
2494
2495 /* clear working vector and reinsert columns to lists
2496 */
2497 for(i = u.row.start[prow], pend = i + plen; i < pend; ++i)
2498 {
2499 j = u.row.idx[i];
2500 work[j] = 0;
2501 temp.s_mark[j] = 0;
2502 init2DR(temp.pivot_col[j], temp.pivot_colNZ[temp.s_cact[j]]);
2503 assert(col.perm[j] < 0);
2504 }
2505 }
2506
2507
2508 /*
2509 * Factorize nucleus.
2510 */
2511 inline void CLUFactorRational::eliminateNucleus(const Rational& threshold)
2512 {
2513 int r, c;
2514 CLUFactorRational::Pring* pivot;
2515
2516 if(stat == SLinSolverRational::SINGULAR)
2517 return;
2518
2519 temp.pivots.mkwtz = -1;
2520
2521 temp.pivots.idx = -1;
2522
2523 temp.pivots.pos = -1;
2524
2525 while(temp.stage < thedim - 1)
2526 {
2527 #ifndef NDEBUG
2528 int i;
2529 // CLUFactorRationalIsConsistent(fac);
2530
2531 for(i = 0; i < thedim; ++i)
2532 if(col.perm[i] < 0)
2533 assert(temp.s_mark[i] == 0);
2534
2535 #endif
2536
2537 if(temp.pivot_rowNZ[1].next != &(temp.pivot_rowNZ[1]))
2538 {
2539 if(timeLimitReached())
2540 return;
2541
2542 /* row singleton available */
2543 eliminateRowSingletons();
2544 }
2545 else if(temp.pivot_colNZ[1].next != &(temp.pivot_colNZ[1]))
2546 {
2547 if(timeLimitReached())
2548 return;
2549
2550 /* column singleton available */
2551 eliminateColSingletons();
2552 }
2553 else
2554 {
2555 initDR(temp.pivots);
2556 selectPivots(threshold);
2557
2558 assert(temp.pivots.next != &temp.pivots &&
2559 "ERROR: no pivot element selected");
2560
2561 for(pivot = temp.pivots.next; pivot != &temp.pivots;
2562 pivot = pivot->next)
2563 {
2564 if(timeLimitReached())
2565 return;
2566
2567 eliminatePivot(pivot->idx, pivot->pos);
2568 }
2569 }
2570
2571 if(temp.pivot_rowNZ->next != temp.pivot_rowNZ ||
2572 temp.pivot_colNZ->next != temp.pivot_colNZ)
2573 {
2574 stat = SLinSolverRational::SINGULAR;
2575 return;
2576 }
2577 }
2578
2579 if(temp.stage < thedim)
2580 {
2581 /* Eliminate remaining element.
2582 * Note, that this must be both, column and row singleton.
2583 */
2584 assert(temp.pivot_rowNZ[1].next != &(temp.pivot_rowNZ[1]) &&
2585 "ERROR: one row must be left");
2586 assert(temp.pivot_colNZ[1].next != &(temp.pivot_colNZ[1]) &&
2587 "ERROR: one col must be left");
2588 r = temp.pivot_rowNZ[1].next->idx;
2589 c = temp.pivot_colNZ[1].next->idx;
2590 u.row.len[r] = 0;
2591 u.col.len[c]--;
2592 setPivot(temp.stage, c, r, u.row.val[u.row.start[r]]);
2593 }
2594 }
2595
2596 /*****************************************************************************/
2597
2598 inline int CLUFactorRational::setupColVals()
2599 {
2600 int i;
2601 int n = thedim;
2602
2603 u.col.val.reDim(u.col.size);
2604
2605 for(i = 0; i < thedim; i++)
2606 u.col.len[i] = 0;
2607
2608 maxabs = 0;
2609
2610 for(i = 0; i < thedim; i++)
2611 {
2612 int k = u.row.start[i];
2613 int* idx = &u.row.idx[k];
2614 Rational* val = &u.row.val[k];
2615 int len = u.row.len[i];
2616
2617 n += len;
2618
2619 while(len-- > 0)
2620 {
2621 assert((*idx >= 0) && (*idx < thedim));
2622
2623 k = u.col.start[*idx] + u.col.len[*idx];
2624
2625 assert((k >= 0) && (k < u.col.size));
2626
2627 u.col.len[*idx]++;
2628
2629 assert(u.col.len[*idx] <= u.col.max[*idx]);
2630
2631 u.col.idx[k] = i;
2632 u.col.val[k] = *val;
2633
2634 if(spxAbs(*val) > maxabs)
2635 maxabs = spxAbs(*val);
2636
2637 idx++;
2638
2639 val++;
2640 }
2641 }
2642
2643 return n;
2644 }
2645
2646 /*****************************************************************************/
2647
2648 #ifdef SOPLEX_WITH_L_ROWS
2649 inline void CLUFactorRational::setupRowVals()
2650 {
2651 int i, j, k, m;
2652 int vecs, mem;
2653 int* l_row;
2654 int* idx;
2655 int* beg;
2656 int* l_ridx;
2657 int* l_rbeg;
2658 int* rorig;
2659 int* rrorig;
2660 int* rperm;
2661 int* rrperm;
2662
2663 vecs = l.firstUpdate;
2664 l_row = l.row;
2665 idx = l.idx;
2666 VectorRational& val = l.val;
2667 int validx = 0;
2668 beg = l.start;
2669 mem = beg[vecs];
2670
2671 if(l.ridx)
2672 spx_free(l.ridx);
2673
2674 if(l.rbeg)
2675 spx_free(l.rbeg);
2676
2677 if(l.rorig)
2678 spx_free(l.rorig);
2679
2680 if(l.rperm)
2681 spx_free(l.rperm);
2682
2683 l.rval.reDim(mem + 1);
2684
2685 spx_alloc(l.ridx, mem);
2686
2687 spx_alloc(l.rbeg, thedim + 1);
2688
2689 spx_alloc(l.rorig, thedim);
2690
2691 spx_alloc(l.rperm, thedim);
2692
2693 l_ridx = l.ridx;
2694
2695 VectorRational& l_rval = l.rval;
2696
2697 l_rbeg = l.rbeg;
2698
2699 rorig = l.rorig;
2700
2701 rrorig = row.orig;
2702
2703 rperm = l.rperm;
2704
2705 rrperm = row.perm;
2706
2707 for(i = thedim; i--; *l_rbeg++ = 0)
2708 {
2709 *rorig++ = *rrorig++;
2710 *rperm++ = *rrperm++;
2711 }
2712
2713 *l_rbeg = 0;
2714
2715 l_rbeg = l.rbeg + 1;
2716
2717 for(i = mem; i--;)
2718 l_rbeg[*idx++]++;
2719
2720 idx = l.idx;
2721
2722 for(m = 0, i = thedim; i--; l_rbeg++)
2723 {
2724 j = *l_rbeg;
2725 *l_rbeg = m;
2726 m += j;
2727 }
2728
2729 assert(m == mem);
2730
2731 l_rbeg = l.rbeg + 1;
2732
2733 for(i = j = 0; i < vecs; ++i)
2734 {
2735 m = l_row[i];
2736 assert(idx == &l.idx[l.start[i]]);
2737
2738 for(; j < beg[i + 1]; j++)
2739 {
2740 k = l_rbeg[*idx++]++;
2741 assert(k < mem);
2742 l_ridx[k] = m;
2743 l_rval[k] = val[validx];
2744 validx++; // was l_rval[k] = *val++; with Rational* val
2745 }
2746 }
2747
2748 assert(l.rbeg[thedim] == mem);
2749
2750 assert(l.rbeg[0] == 0);
2751 }
2752
2753 #endif
2754
2755 /*****************************************************************************/
2756
2757 inline void CLUFactorRational::factor(
2758 const SVectorRational** vec, ///< Array of column vector pointers
2759 const Rational& threshold ///< pivoting threshold
2760 )
2761 {
2762 SPxOut::debug(this, "CLUFactorRational::factor()\n");
2763
2764 factorTime->start();
2765
2766 stat = SLinSolverRational::OK;
2767
2768 l.start[0] = 0;
2769 l.firstUpdate = 0;
2770 l.firstUnused = 0;
2771
2772 temp.init(thedim);
2773 initPerm();
2774
2775 initFactorMatrix(vec);
2776
2777 if(stat)
2778 goto TERMINATE;
2779
2780 if(timeLimitReached())
2781 goto TERMINATE;
2782
2783 // initMaxabs = initMaxabs;
2784
2785 colSingletons();
2786
2787 if(stat != SLinSolverRational::OK)
2788 goto TERMINATE;
2789
2790 if(timeLimitReached())
2791 goto TERMINATE;
2792
2793 rowSingletons();
2794
2795 if(stat != SLinSolverRational::OK)
2796 goto TERMINATE;
2797
2798 if(temp.stage < thedim)
2799 {
2800 if(timeLimitReached())
2801 goto TERMINATE;
2802
2803 initFactorRings();
2804 eliminateNucleus(threshold);
2805 freeFactorRings();
2806 }
2807
2808 TERMINATE:
2809
2810 l.firstUpdate = l.firstUnused;
2811
2812 if(stat == SLinSolverRational::OK)
2813 {
2814 #ifdef SOPLEX_WITH_L_ROWS
2815 setupRowVals();
2816 #endif
2817 nzCnt = setupColVals();
2818 }
2819
2820 factorTime->stop();
2821
2822 factorCount++;
2823 }
2824
2825 inline void CLUFactorRational::dump() const
2826 {
2827 int i, j, k;
2828
2829 // Dump regardless of the verbosity level if this method is called;
2830
2831 /* Dump U:
2832 */
2833
2834 for(i = 0; i < thedim; ++i)
2835 {
2836 if(row.perm[i] >= 0)
2837 std::cout << "DCLUFA01 diag[" << i << "]: [" << col.orig[row.perm[i]]
2838 << "] = " << diag[i] << std::endl;
2839
2840 for(j = 0; j < u.row.len[i]; ++j)
2841 std::cout << "DCLUFA02 u[" << i << "]: ["
2842 << u.row.idx[u.row.start[i] + j] << "] = "
2843 << u.row.val[u.row.start[i] + j] << std::endl;
2844 }
2845
2846 /* Dump L:
2847 */
2848 for(i = 0; i < thedim; ++i)
2849 {
2850 for(j = 0; j < l.firstUnused; ++j)
2851 if(col.orig[row.perm[l.row[j]]] == i)
2852 {
2853 std::cout << "DCLUFA03 l[" << i << "]" << std::endl;
2854
2855 for(k = l.start[j]; k < l.start[j + 1]; ++k)
2856 std::cout << "DCLUFA04 l[" << k - l.start[j]
2857 << "]: [" << l.idx[k]
2858 << "] = " << l.val[k] << std::endl;
2859
2860 break;
2861 }
2862 }
2863
2864 return;
2865 }
2866
2867 /*****************************************************************************/
2868 /*
2869 * Ensure that row memory is at least size.
2870 */
2871 inline void CLUFactorRational::minRowMem(int size)
2872 {
2873 if(u.row.val.dim() < size)
2874 {
2875 u.row.val.reDim(size);
2876 spx_realloc(u.row.idx, size);
2877 }
2878
2879 assert(u.row.val.dim() >= size);
2880 }
2881
2882 /*****************************************************************************/
2883 /*
2884 * Ensure that column memory is at least size.
2885 */
2886 inline void CLUFactorRational::minColMem(int size)
2887 {
2888
2889 if(u.col.size < size)
2890 {
2891 u.col.size = size;
2892 spx_realloc(u.col.idx, size);
2893 }
2894 }
2895
2896 inline void CLUFactorRational::forestMinColMem(int size)
2897 {
2898
2899 if(u.col.size < size)
2900 {
2901 u.col.size = size;
2902 spx_realloc(u.col.idx, size);
2903 u.col.val.reDim(size);
2904 }
2905 }
2906
2907 inline void CLUFactorRational::minLMem(int size)
2908 {
2909
2910 if(size > l.val.dim())
2911 {
2912 int newsize = int(0.2 * l.val.dim() + size);
2913 l.val.reDim(newsize);
2914 spx_realloc(l.idx, l.val.dim());
2915 }
2916 }
2917
2918
2919 inline int CLUFactorRational::makeLvec(int p_len, int p_row)
2920 {
2921
2922 if(l.firstUnused >= l.startSize)
2923 {
2924 l.startSize += 100;
2925 spx_realloc(l.start, l.startSize);
2926 }
2927
2928 int* p_lrow = l.row;
2929
2930 int* p_lbeg = l.start;
2931 int first = p_lbeg[l.firstUnused];
2932
2933 assert(p_len > 0 && "ERROR: no empty columns allowed in L vectors");
2934
2935 minLMem(first + p_len);
2936 p_lrow[l.firstUnused] = p_row;
2937 l.start[++(l.firstUnused)] = first + p_len;
2938
2939 assert(l.start[l.firstUnused] <= l.val.dim());
2940 assert(l.firstUnused <= l.startSize);
2941 return first;
2942 }
2943
2944
2945 /*****************************************************************************/
2946
2947 inline bool CLUFactorRational::isConsistent() const
2948 {
2949 #ifdef ENABLE_CONSISTENCY_CHECKS
2950 int i, j, k, ll;
2951 Dring* ring;
2952 CLUFactorRational::Pring* pring;
2953
2954 /* Consistency only relevant for real factorizations
2955 */
2956
2957 if(stat)
2958 return true;
2959
2960 /* Test column ring list consistency.
2961 */
2962 i = 0;
2963
2964 for(ring = u.col.list.next; ring != &(u.col.list); ring = ring->next)
2965 {
2966 assert(ring->idx >= 0);
2967 assert(ring->idx < thedim);
2968 assert(ring->prev->next == ring);
2969
2970 if(ring != u.col.list.next)
2971 {
2972 assert(u.col.start[ring->prev->idx] + u.col.max[ring->prev->idx]
2973 == u.col.start[ring->idx]);
2974 }
2975
2976 ++i;
2977 }
2978
2979 assert(i == thedim);
2980
2981 assert(u.col.start[ring->prev->idx] + u.col.max[ring->prev->idx]
2982 == u.col.used);
2983
2984
2985 /* Test row ring list consistency.
2986 */
2987 i = 0;
2988
2989 for(ring = u.row.list.next; ring != &(u.row.list); ring = ring->next)
2990 {
2991 assert(ring->idx >= 0);
2992 assert(ring->idx < thedim);
2993 assert(ring->prev->next == ring);
2994 assert(u.row.start[ring->prev->idx] + u.row.max[ring->prev->idx]
2995 == u.row.start[ring->idx]);
2996 ++i;
2997 }
2998
2999 assert(i == thedim);
3000
3001 assert(u.row.start[ring->prev->idx] + u.row.max[ring->prev->idx]
3002 == u.row.used);
3003
3004
3005 /* Test consistency of individual svectors.
3006 */
3007
3008 for(i = 0; i < thedim; ++i)
3009 {
3010 assert(u.row.max[i] >= u.row.len[i]);
3011 assert(u.col.max[i] >= u.col.len[i]);
3012 }
3013
3014
3015 /* Test consistency of column file to row file of U
3016 */
3017 for(i = 0; i < thedim; ++i)
3018 {
3019 for(j = u.row.start[i] + u.row.len[i] - 1; j >= u.row.start[i]; j--)
3020 {
3021 k = u.row.idx[j];
3022
3023 for(ll = u.col.start[k] + u.col.len[k] - 1; ll >= u.col.start[k]; ll--)
3024 {
3025 if(u.col.idx[ll] == i)
3026 break;
3027 }
3028
3029 assert(u.col.idx[ll] == i);
3030
3031 if(row.perm[i] < 0)
3032 {
3033 assert(col.perm[k] < 0);
3034 }
3035 else
3036 {
3037 assert(col.perm[k] < 0 || col.perm[k] > row.perm[i]);
3038 }
3039 }
3040 }
3041
3042 /* Test consistency of row file to column file of U
3043 */
3044 for(i = 0; i < thedim; ++i)
3045 {
3046 for(j = u.col.start[i] + u.col.len[i] - 1; j >= u.col.start[i]; j--)
3047 {
3048 k = u.col.idx[j];
3049
3050 for(ll = u.row.start[k] + u.row.len[k] - 1; ll >= u.row.start[k]; ll--)
3051 {
3052 if(u.row.idx[ll] == i)
3053 break;
3054 }
3055
3056 assert(u.row.idx[ll] == i);
3057
3058 assert(col.perm[i] < 0 || row.perm[k] < col.perm[i]);
3059 }
3060 }
3061
3062 /* Test consistency of nonzero count lists
3063 */
3064 if(temp.pivot_colNZ && temp.pivot_rowNZ)
3065 {
3066 for(i = 0; i < thedim - temp.stage; ++i)
3067 {
3068 for(pring = temp.pivot_rowNZ[i].next; pring != &(temp.pivot_rowNZ[i]); pring = pring->next)
3069 {
3070 assert(row.perm[pring->idx] < 0);
3071 }
3072
3073 for(pring = temp.pivot_colNZ[i].next; pring != &(temp.pivot_colNZ[i]); pring = pring->next)
3074 {
3075 assert(col.perm[pring->idx] < 0);
3076 }
3077 }
3078 }
3079
3080 #endif // CONSISTENCY_CHECKS
3081
3082 return true;
3083 }
3084
3085 inline void CLUFactorRational::solveUright(Rational* wrk, Rational* vec)
3086 {
3087
3088 for(int i = thedim - 1; i >= 0; i--)
3089 {
3090 int r = row.orig[i];
3091 int c = col.orig[i];
3092 Rational x = wrk[c] = diag[r] * vec[r];
3093
3094 vec[r] = 0;
3095
3096 if(x != 0)
3097 //if (isNotZero(x))
3098 {
3099 if(timeLimitReached())
3100 return;
3101
3102 for(int j = u.col.start[c]; j < u.col.start[c] + u.col.len[c]; j++)
3103 vec[u.col.idx[j]] -= x * u.col.val[j];
3104 }
3105 }
3106 }
3107
3108 inline int CLUFactorRational::solveUrightEps(Rational* vec, int* nonz, Rational* rhs)
3109 {
3110 int i, j, r, c, n;
3111 int* rorig, *corig;
3112 int* cidx, *clen, *cbeg;
3113 Rational x;
3114
3115 int* idx;
3116 Rational* val;
3117
3118 rorig = row.orig;
3119 corig = col.orig;
3120
3121 cidx = u.col.idx;
3122 VectorRational& cval = u.col.val;
3123 clen = u.col.len;
3124 cbeg = u.col.start;
3125
3126 n = 0;
3127
3128 for(i = thedim - 1; i >= 0; --i)
3129 {
3130 r = rorig[i];
3131 x = diag[r] * rhs[r];
3132
3133 if(x != 0)
3134 {
3135 c = corig[i];
3136 vec[c] = x;
3137 nonz[n++] = c;
3138 val = &cval[cbeg[c]];
3139 idx = &cidx[cbeg[c]];
3140 j = clen[c];
3141
3142 while(j-- > 0)
3143 rhs[*idx++] -= x * (*val++);
3144 }
3145 }
3146
3147 return n;
3148 }
3149
3150 inline void CLUFactorRational::solveUright2(Rational* p_work1, Rational* vec1, Rational* p_work2,
3151 Rational* vec2)
3152 {
3153 int i, j, r, c;
3154 int* rorig, *corig;
3155 int* cidx, *clen, *cbeg;
3156 Rational x1, x2;
3157
3158 int* idx;
3159 Rational* val;
3160
3161 rorig = row.orig;
3162 corig = col.orig;
3163
3164 cidx = u.col.idx;
3165 VectorRational& cval = u.col.val;
3166 clen = u.col.len;
3167 cbeg = u.col.start;
3168
3169 for(i = thedim - 1; i >= 0; --i)
3170 {
3171 r = rorig[i];
3172 c = corig[i];
3173 p_work1[c] = x1 = diag[r] * vec1[r];
3174 p_work2[c] = x2 = diag[r] * vec2[r];
3175 vec1[r] = vec2[r] = 0;
3176
3177 if(x1 != 0 && x2 != 0)
3178 {
3179 val = &cval[cbeg[c]];
3180 idx = &cidx[cbeg[c]];
3181 j = clen[c];
3182
3183 while(j-- > 0)
3184 {
3185 vec1[*idx] -= x1 * (*val);
3186 vec2[*idx++] -= x2 * (*val++);
3187 }
3188 }
3189 else if(x1 != 0)
3190 {
3191 val = &cval[cbeg[c]];
3192 idx = &cidx[cbeg[c]];
3193 j = clen[c];
3194
3195 while(j-- > 0)
3196 vec1[*idx++] -= x1 * (*val++);
3197 }
3198 else if(x2 != 0)
3199 {
3200 val = &cval[cbeg[c]];
3201 idx = &cidx[cbeg[c]];
3202 j = clen[c];
3203
3204 while(j-- > 0)
3205 vec2[*idx++] -= x2 * (*val++);
3206 }
3207 }
3208 }
3209
3210 inline int CLUFactorRational::solveUright2eps(Rational* p_work1, Rational* vec1, Rational* p_work2,
3211 Rational* vec2, int* nonz)
3212 {
3213 int i, j, r, c, n;
3214 int* rorig, *corig;
3215 int* cidx, *clen, *cbeg;
3216 bool notzero1, notzero2;
3217 Rational x1, x2;
3218
3219 int* idx;
3220 Rational* val;
3221
3222 rorig = row.orig;
3223 corig = col.orig;
3224
3225 cidx = u.col.idx;
3226 VectorRational& cval = u.col.val;
3227 clen = u.col.len;
3228 cbeg = u.col.start;
3229
3230 n = 0;
3231
3232 for(i = thedim - 1; i >= 0; --i)
3233 {
3234 c = corig[i];
3235 r = rorig[i];
3236 p_work1[c] = x1 = diag[r] * vec1[r];
3237 p_work2[c] = x2 = diag[r] * vec2[r];
3238 vec1[r] = vec2[r] = 0;
3239 notzero1 = x1 != 0 ? 1 : 0;
3240 notzero2 = x2 != 0 ? 1 : 0;
3241
3242 if(notzero1 && notzero2)
3243 {
3244 *nonz++ = c;
3245 n++;
3246 val = &cval[cbeg[c]];
3247 idx = &cidx[cbeg[c]];
3248 j = clen[c];
3249
3250 while(j-- > 0)
3251 {
3252 vec1[*idx] -= x1 * (*val);
3253 vec2[*idx++] -= x2 * (*val++);
3254 }
3255 }
3256 else if(notzero1)
3257 {
3258 p_work2[c] = 0;
3259 *nonz++ = c;
3260 n++;
3261 val = &cval[cbeg[c]];
3262 idx = &cidx[cbeg[c]];
3263 j = clen[c];
3264
3265 while(j-- > 0)
3266 vec1[*idx++] -= x1 * (*val++);
3267 }
3268 else if(notzero2)
3269 {
3270 p_work1[c] = 0;
3271 val = &cval[cbeg[c]];
3272 idx = &cidx[cbeg[c]];
3273 j = clen[c];
3274
3275 while(j-- > 0)
3276 vec2[*idx++] -= x2 * (*val++);
3277 }
3278 else
3279 {
3280 p_work1[c] = 0;
3281 p_work2[c] = 0;
3282 }
3283 }
3284
3285 return n;
3286 }
3287
3288 inline void CLUFactorRational::solveLright(Rational* vec)
3289 {
3290 int i, j, k;
3291 int end;
3292 Rational x;
3293 Rational* val;
3294 int* lrow, *lidx, *idx;
3295 int* lbeg;
3296
3297 VectorRational& lval = l.val;
3298 lidx = l.idx;
3299 lrow = l.row;
3300 lbeg = l.start;
3301
3302 end = l.firstUpdate;
3303
3304 for(i = 0; i < end; ++i)
3305 {
3306 if((x = vec[lrow[i]]) != 0)
3307 {
3308 if(timeLimitReached())
3309 return;
3310
3311 SPxOut::debug(this, "y{}={}\n", lrow[i], vec[lrow[i]]);
3312
3313 k = lbeg[i];
3314 idx = &(lidx[k]);
3315 val = &(lval[k]);
3316
3317 for(j = lbeg[i + 1]; j > k; --j)
3318 {
3319 SPxOut::debug(this, " -> y{} -= {} * {} = {} -> {}\n", *idx, x, *val,
3320 x * (*val), vec[*idx] - x * (*val));
3321 vec[*idx++] -= x * (*val++);
3322 }
3323 }
3324 }
3325
3326 if(l.updateType) /* Forest-Tomlin Updates */
3327 {
3328 SPxOut::debug(this, "performing FT updates...\n");
3329
3330 end = l.firstUnused;
3331
3332 for(; i < end; ++i)
3333 {
3334 x = 0;
3335 k = lbeg[i];
3336 idx = &(lidx[k]);
3337 val = &(lval[k]);
3338
3339 for(j = lbeg[i + 1]; j > k; --j)
3340 x += vec[*idx++] * (*val++);
3341
3342 vec[lrow[i]] -= x;
3343
3344 SPxOut::debug(this, "y{}={}\n", lrow[i], vec[lrow[i]]);
3345 }
3346
3347 SPxOut::debug(this, "finished FT updates.\n");
3348 }
3349 }
3350
3351 inline void CLUFactorRational::solveLright2(Rational* vec1, Rational* vec2)
3352 {
3353 int i, j, k;
3354 int end;
3355 Rational x2;
3356 Rational x1;
3357 Rational* val;
3358 int* lrow, *lidx, *idx;
3359 int* lbeg;
3360
3361 VectorRational& lval = l.val;
3362 lidx = l.idx;
3363 lrow = l.row;
3364 lbeg = l.start;
3365
3366 end = l.firstUpdate;
3367
3368 for(i = 0; i < end; ++i)
3369 {
3370 x1 = vec1[lrow[i]];
3371 x2 = vec2[lrow[i]];
3372
3373 if(x1 != 0 && x2 != 0)
3374 {
3375 k = lbeg[i];
3376 idx = &(lidx[k]);
3377 val = &(lval[k]);
3378
3379 for(j = lbeg[i + 1]; j > k; --j)
3380 {
3381 vec1[*idx] -= x1 * (*val);
3382 vec2[*idx++] -= x2 * (*val++);
3383 }
3384 }
3385 else if(x1 != 0)
3386 {
3387 k = lbeg[i];
3388 idx = &(lidx[k]);
3389 val = &(lval[k]);
3390
3391 for(j = lbeg[i + 1]; j > k; --j)
3392 vec1[*idx++] -= x1 * (*val++);
3393 }
3394 else if(x2 != 0)
3395 {
3396 k = lbeg[i];
3397 idx = &(lidx[k]);
3398 val = &(lval[k]);
3399
3400 for(j = lbeg[i + 1]; j > k; --j)
3401 vec2[*idx++] -= x2 * (*val++);
3402 }
3403 }
3404
3405 if(l.updateType) /* Forest-Tomlin Updates */
3406 {
3407 end = l.firstUnused;
3408
3409 for(; i < end; ++i)
3410 {
3411 x1 = 0;
3412 x2 = 0;
3413 k = lbeg[i];
3414 idx = &(lidx[k]);
3415 val = &(lval[k]);
3416
3417 for(j = lbeg[i + 1]; j > k; --j)
3418 {
3419 x1 += vec1[*idx] * (*val);
3420 x2 += vec2[*idx++] * (*val++);
3421 }
3422
3423 vec1[lrow[i]] -= x1;
3424
3425 vec2[lrow[i]] -= x2;
3426 }
3427 }
3428 }
3429
3430 inline void CLUFactorRational::solveUpdateRight(Rational* vec)
3431 {
3432 int i, j, k;
3433 int end;
3434 Rational x;
3435 Rational* val;
3436 int* lrow, *lidx, *idx;
3437 int* lbeg;
3438
3439 assert(!l.updateType); /* no Forest-Tomlin Updates */
3440
3441 VectorRational& lval = l.val;
3442 lidx = l.idx;
3443 lrow = l.row;
3444 lbeg = l.start;
3445
3446 end = l.firstUnused;
3447
3448 for(i = l.firstUpdate; i < end; ++i)
3449 {
3450 if((x = vec[lrow[i]]) != 0)
3451 {
3452 k = lbeg[i];
3453 idx = &(lidx[k]);
3454 val = &(lval[k]);
3455
3456 for(j = lbeg[i + 1]; j > k; --j)
3457 vec[*idx++] -= x * (*val++);
3458 }
3459 }
3460 }
3461
3462 inline void CLUFactorRational::solveUpdateRight2(Rational* vec1, Rational* vec2)
3463 {
3464 int i, j, k;
3465 int end;
3466 Rational x1, x2;
3467 int* lrow, *lidx;
3468 int* lbeg;
3469
3470 int* idx;
3471 Rational* val;
3472
3473 assert(!l.updateType); /* no Forest-Tomlin Updates */
3474
3475 VectorRational& lval = l.val;
3476 lidx = l.idx;
3477 lrow = l.row;
3478 lbeg = l.start;
3479
3480 end = l.firstUnused;
3481
3482 for(i = l.firstUpdate; i < end; ++i)
3483 {
3484 x1 = vec1[lrow[i]];
3485 x2 = vec2[lrow[i]];
3486
3487 if(x1 != 0 && x2 != 0)
3488 {
3489 k = lbeg[i];
3490 idx = &(lidx[k]);
3491 val = &(lval[k]);
3492
3493 for(j = lbeg[i + 1]; j > k; --j)
3494 {
3495 vec1[*idx] -= x1 * (*val);
3496 vec2[*idx++] -= x2 * (*val++);
3497 }
3498 }
3499 else if(x1 != 0)
3500 {
3501 k = lbeg[i];
3502 idx = &(lidx[k]);
3503 val = &(lval[k]);
3504
3505 for(j = lbeg[i + 1]; j > k; --j)
3506 vec1[*idx++] -= x1 * (*val++);
3507 }
3508 else if(x2 != 0)
3509 {
3510 k = lbeg[i];
3511 idx = &(lidx[k]);
3512 val = &(lval[k]);
3513
3514 for(j = lbeg[i + 1]; j > k; --j)
3515 vec2[*idx++] -= x2 * (*val++);
3516 }
3517 }
3518 }
3519
3520 inline int CLUFactorRational::solveRight4update(Rational* vec, int* nonz,
3521 Rational* rhs, Rational* forest, int* forestNum, int* forestIdx)
3522 {
3523 solveLright(rhs);
3524
3525 if(forest)
3526 {
3527 int n = 0;
3528
3529 for(int i = 0; i < thedim; i++)
3530 {
3531 forestIdx[n] = i;
3532 forest[i] = rhs[i];
3533 n += rhs[i] != 0 ? 1 : 0;
3534 }
3535
3536 *forestNum = n;
3537 }
3538
3539 if(!l.updateType) /* no Forest-Tomlin Updates */
3540 {
3541 solveUright(vec, rhs);
3542 solveUpdateRight(vec);
3543 return 0;
3544 }
3545 else
3546 return solveUrightEps(vec, nonz, rhs);
3547 }
3548
3549 inline void CLUFactorRational::solveRight(Rational* vec, Rational* rhs)
3550 {
3551 solveLright(rhs);
3552 solveUright(vec, rhs);
3553
3554 if(!l.updateType) /* no Forest-Tomlin Updates */
3555 solveUpdateRight(vec);
3556 }
3557
3558 inline int CLUFactorRational::solveRight2update(Rational* vec1,
3559 Rational* vec2,
3560 Rational* rhs1,
3561 Rational* rhs2,
3562 int* nonz,
3563 Rational* forest,
3564 int* forestNum,
3565 int* forestIdx)
3566 {
3567 solveLright2(rhs1, rhs2);
3568
3569 if(forest)
3570 {
3571 int n = 0;
3572
3573 for(int i = 0; i < thedim; i++)
3574 {
3575 forestIdx[n] = i;
3576 forest[i] = rhs1[i];
3577 n += rhs1[i] != 0 ? 1 : 0;
3578 }
3579
3580 *forestNum = n;
3581 }
3582
3583 if(!l.updateType) /* no Forest-Tomlin Updates */
3584 {
3585 solveUright2(vec1, rhs1, vec2, rhs2);
3586 solveUpdateRight2(vec1, vec2);
3587 return 0;
3588 }
3589 else
3590 return solveUright2eps(vec1, rhs1, vec2, rhs2, nonz);
3591 }
3592
3593 inline void CLUFactorRational::solveRight2(
3594 Rational* vec1,
3595 Rational* vec2,
3596 Rational* rhs1,
3597 Rational* rhs2)
3598 {
3599 solveLright2(rhs1, rhs2);
3600
3601 if(l.updateType) /* Forest-Tomlin Updates */
3602 solveUright2(vec1, rhs1, vec2, rhs2);
3603 else
3604 {
3605 solveUright2(vec1, rhs1, vec2, rhs2);
3606 solveUpdateRight2(vec1, vec2);
3607 }
3608 }
3609
3610 /*****************************************************************************/
3611 inline void CLUFactorRational::solveUleft(Rational* p_work, Rational* vec)
3612 {
3613 for(int i = 0; i < thedim; ++i)
3614 {
3615 int c = col.orig[i];
3616 int r = row.orig[i];
3617
3618 assert(c >= 0); // Inna/Tobi: otherwise, vec[c] would be strange...
3619 assert(r >= 0); // Inna/Tobi: otherwise, diag[r] would be strange...
3620
3621 Rational x = vec[c];
3622
3623 vec[c] = 0;
3624
3625 if(x != 0)
3626 {
3627 if(timeLimitReached())
3628 return;
3629
3630 x *= diag[r];
3631 p_work[r] = x;
3632
3633 int end = u.row.start[r] + u.row.len[r];
3634
3635 for(int m = u.row.start[r]; m < end; m++)
3636 {
3637 vec[u.row.idx[m]] -= x * u.row.val[m];
3638 }
3639 }
3640 }
3641 }
3642
3643
3644
3645 inline void CLUFactorRational::solveUleft2(Rational* p_work1, Rational* vec1, Rational* p_work2,
3646 Rational* vec2)
3647 {
3648 Rational x1;
3649 Rational x2;
3650 int i, k, r, c;
3651 int* rorig, *corig;
3652 int* ridx, *rlen, *rbeg, *idx;
3653 Rational* val;
3654
3655 rorig = row.orig;
3656 corig = col.orig;
3657
3658 ridx = u.row.idx;
3659 VectorRational& rval = u.row.val;
3660 rlen = u.row.len;
3661 rbeg = u.row.start;
3662
3663 for(i = 0; i < thedim; ++i)
3664 {
3665 c = corig[i];
3666 r = rorig[i];
3667 x1 = vec1[c];
3668 x2 = vec2[c];
3669
3670 if((x1 != 0) && (x2 != 0))
3671 {
3672 x1 *= diag[r];
3673 x2 *= diag[r];
3674 p_work1[r] = x1;
3675 p_work2[r] = x2;
3676 k = rbeg[r];
3677 idx = &ridx[k];
3678 val = &rval[k];
3679
3680 for(int m = rlen[r]; m != 0; --m)
3681 {
3682 vec1[*idx] -= x1 * (*val);
3683 vec2[*idx] -= x2 * (*val++);
3684 idx++;
3685 }
3686 }
3687 else if(x1 != 0)
3688 {
3689 x1 *= diag[r];
3690 p_work1[r] = x1;
3691 k = rbeg[r];
3692 idx = &ridx[k];
3693 val = &rval[k];
3694
3695 for(int m = rlen[r]; m != 0; --m)
3696 vec1[*idx++] -= x1 * (*val++);
3697 }
3698 else if(x2 != 0)
3699 {
3700 x2 *= diag[r];
3701 p_work2[r] = x2;
3702 k = rbeg[r];
3703 idx = &ridx[k];
3704 val = &rval[k];
3705
3706 for(int m = rlen[r]; m != 0; --m)
3707 vec2[*idx++] -= x2 * (*val++);
3708 }
3709 }
3710 }
3711
3712 inline int CLUFactorRational::solveLleft2forest(Rational* vec1, int* /* nonz */, Rational* vec2)
3713 {
3714 int i;
3715 int j;
3716 int k;
3717 int end;
3718 Rational x1, x2;
3719 Rational* val;
3720 int* lidx, *idx, *lrow;
3721 int* lbeg;
3722
3723 VectorRational& lval = l.val;
3724 lidx = l.idx;
3725 lrow = l.row;
3726 lbeg = l.start;
3727
3728 end = l.firstUpdate;
3729
3730 for(i = l.firstUnused - 1; i >= end; --i)
3731 {
3732 j = lrow[i];
3733 x1 = vec1[j];
3734 x2 = vec2[j];
3735
3736 if(x1 != 0)
3737 {
3738 if(x2 != 0)
3739 {
3740 k = lbeg[i];
3741 val = &lval[k];
3742 idx = &lidx[k];
3743
3744 for(j = lbeg[i + 1]; j > k; --j)
3745 {
3746 vec1[*idx] -= x1 * (*val);
3747 vec2[*idx++] -= x2 * (*val++);
3748 }
3749 }
3750 else
3751 {
3752 k = lbeg[i];
3753 val = &lval[k];
3754 idx = &lidx[k];
3755
3756 for(j = lbeg[i + 1]; j > k; --j)
3757 vec1[*idx++] -= x1 * (*val++);
3758 }
3759 }
3760 else if(x2 != 0)
3761 {
3762 k = lbeg[i];
3763 val = &lval[k];
3764 idx = &lidx[k];
3765
3766 for(j = lbeg[i + 1]; j > k; --j)
3767 vec2[*idx++] -= x2 * (*val++);
3768 }
3769 }
3770
3771 return 0;
3772 }
3773
3774 inline void CLUFactorRational::solveLleft2(Rational* vec1, int* /* nonz */, Rational* vec2)
3775 {
3776 int i, j, k, r;
3777 int x1not0, x2not0;
3778 Rational x1, x2;
3779
3780 Rational* val;
3781 int* ridx, *idx;
3782 int* rbeg;
3783 int* rorig;
3784
3785 ridx = l.ridx;
3786 VectorRational& rval = l.rval;
3787 rbeg = l.rbeg;
3788 rorig = l.rorig;
3789
3790 #ifndef SOPLEX_WITH_L_ROWS
3791 VectorRational& lval = l.val;
3792 int* lidx = l.idx;
3793 int* lrow = l.row;
3794 int* lbeg = l.start;
3795
3796 i = l.firstUpdate - 1;
3797
3798 for(; i >= 0; --i)
3799 {
3800 k = lbeg[i];
3801 val = &lval[k];
3802 idx = &lidx[k];
3803 x1 = 0;
3804 x2 = 0;
3805
3806 for(j = lbeg[i + 1]; j > k; --j)
3807 {
3808 x1 += vec1[*idx] * (*val);
3809 x2 += vec2[*idx++] * (*val++);
3810 }
3811
3812 vec1[lrow[i]] -= x1;
3813
3814 vec2[lrow[i]] -= x2;
3815 }
3816
3817 #else
3818
3819 for(i = thedim; i--;)
3820 {
3821 r = rorig[i];
3822 x1 = vec1[r];
3823 x2 = vec2[r];
3824 x1not0 = (x1 != 0);
3825 x2not0 = (x2 != 0);
3826
3827 if(x1not0 && x2not0)
3828 {
3829 k = rbeg[r];
3830 j = rbeg[r + 1] - k;
3831 val = &rval[k];
3832 idx = &ridx[k];
3833
3834 while(j-- > 0)
3835 {
3836 assert(row.perm[*idx] < i);
3837 vec1[*idx] -= x1 * *val;
3838 vec2[*idx++] -= x2 * *val++;
3839 }
3840 }
3841 else if(x1not0)
3842 {
3843 k = rbeg[r];
3844 j = rbeg[r + 1] - k;
3845 val = &rval[k];
3846 idx = &ridx[k];
3847
3848 while(j-- > 0)
3849 {
3850 assert(row.perm[*idx] < i);
3851 vec1[*idx++] -= x1 * *val++;
3852 }
3853 }
3854 else if(x2not0)
3855 {
3856 k = rbeg[r];
3857 j = rbeg[r + 1] - k;
3858 val = &rval[k];
3859 idx = &ridx[k];
3860
3861 while(j-- > 0)
3862 {
3863 assert(row.perm[*idx] < i);
3864 vec2[*idx++] -= x2 * *val++;
3865 }
3866 }
3867 }
3868
3869 #endif
3870 }
3871
3872 inline int CLUFactorRational::solveLleftForest(Rational* vec, int* /* nonz */)
3873 {
3874 int i, j, k, end;
3875 Rational x;
3876 Rational* val;
3877 int* idx, *lidx, *lrow, *lbeg;
3878
3879 VectorRational& lval = l.val;
3880 lidx = l.idx;
3881 lrow = l.row;
3882 lbeg = l.start;
3883
3884 end = l.firstUpdate;
3885
3886 for(i = l.firstUnused - 1; i >= end; --i)
3887 {
3888 if((x = vec[lrow[i]]) != 0)
3889 {
3890 if(timeLimitReached())
3891 return 0;
3892
3893 k = lbeg[i];
3894 val = &lval[k];
3895 idx = &lidx[k];
3896
3897 for(j = lbeg[i + 1]; j > k; --j)
3898 vec[*idx++] -= x * (*val++);
3899 }
3900 }
3901
3902 return 0;
3903 }
3904
3905 inline void CLUFactorRational::solveLleft(Rational* vec)
3906 {
3907
3908 #ifndef SOPLEX_WITH_L_ROWS
3909 int* idx;
3910 Rational* val;
3911 VectorRational& lval = l.val;
3912 int* lidx = l.idx;
3913 int* lrow = l.row;
3914 int* lbeg = l.start;
3915
3916 for(int i = l.firstUpdate - 1; i >= 0; --i)
3917 {
3918 if(timeLimitReached())
3919 return;
3920
3921 int k = lbeg[i];
3922 val = &lval[k];
3923 idx = &lidx[k];
3924 x = 0;
3925
3926 for(int j = lbeg[i + 1]; j > k; --j)
3927 x += vec[*idx++] * (*val++);
3928
3929 vec[lrow[i]] -= x;
3930 }
3931
3932 #else
3933
3934 for(int i = thedim - 1; i >= 0; --i)
3935 {
3936 int r = l.rorig[i];
3937 Rational x = vec[r];
3938
3939 if(x != 0)
3940 {
3941 if(timeLimitReached())
3942 return;
3943
3944 for(int k = l.rbeg[r]; k < l.rbeg[r + 1]; k++)
3945 {
3946 int j = l.ridx[k];
3947
3948 assert(l.rperm[j] < i);
3949
3950 vec[j] -= x * l.rval[k];
3951 }
3952 }
3953 }
3954
3955 #endif // SOPLEX_WITH_L_ROWS
3956 }
3957
3958 inline int CLUFactorRational::solveLleftEps(Rational* vec, int* nonz)
3959 {
3960 int i, j, k, n;
3961 int r;
3962 Rational x;
3963 Rational* val;
3964 int* ridx, *idx;
3965 int* rbeg;
3966 int* rorig;
3967
3968 ridx = l.ridx;
3969 VectorRational& rval = l.rval;
3970 rbeg = l.rbeg;
3971 rorig = l.rorig;
3972 n = 0;
3973 #ifndef SOPLEX_WITH_L_ROWS
3974 VectorRational& lval = l.val;
3975 int* lidx = l.idx;
3976 int* lrow = l.row;
3977 int* lbeg = l.start;
3978
3979 for(i = l.firstUpdate - 1; i >= 0; --i)
3980 {
3981 k = lbeg[i];
3982 val = &lval[k];
3983 idx = &lidx[k];
3984 x = 0;
3985
3986 for(j = lbeg[i + 1]; j > k; --j)
3987 x += vec[*idx++] * (*val++);
3988
3989 vec[lrow[i]] -= x;
3990 }
3991
3992 #else
3993
3994 for(i = thedim; i--;)
3995 {
3996 r = rorig[i];
3997 x = vec[r];
3998
3999 if(x != 0)
4000 {
4001 *nonz++ = r;
4002 n++;
4003 k = rbeg[r];
4004 j = rbeg[r + 1] - k;
4005 val = &rval[k];
4006 idx = &ridx[k];
4007
4008 while(j-- > 0)
4009 {
4010 assert(row.perm[*idx] < i);
4011 vec[*idx++] -= x * *val++;
4012 }
4013 }
4014 else
4015 vec[r] = 0;
4016 }
4017
4018 #endif
4019
4020 return n;
4021 }
4022
4023 inline void CLUFactorRational::solveUpdateLeft(Rational* vec)
4024 {
4025 int i, j, k, end;
4026 Rational x;
4027 Rational* val;
4028 int* lrow, *lidx, *idx;
4029 int* lbeg;
4030
4031 VectorRational& lval = l.val;
4032 lidx = l.idx;
4033 lrow = l.row;
4034 lbeg = l.start;
4035
4036 assert(!l.updateType); /* Forest-Tomlin Updates */
4037
4038 end = l.firstUpdate;
4039
4040 for(i = l.firstUnused - 1; i >= end; --i)
4041 {
4042 k = lbeg[i];
4043 val = &lval[k];
4044 idx = &lidx[k];
4045 x = 0;
4046
4047 for(j = lbeg[i + 1]; j > k; --j)
4048 x += vec[*idx++] * (*val++);
4049
4050 vec[lrow[i]] -= x;
4051 }
4052 }
4053
4054 inline void CLUFactorRational::solveUpdateLeft2(Rational* vec1, Rational* vec2)
4055 {
4056 int i, j, k, end;
4057 Rational x1, x2;
4058 Rational* val;
4059 int* lrow, *lidx, *idx;
4060 int* lbeg;
4061
4062 VectorRational& lval = l.val;
4063 lidx = l.idx;
4064 lrow = l.row;
4065 lbeg = l.start;
4066
4067 assert(!l.updateType); /* Forest-Tomlin Updates */
4068
4069 end = l.firstUpdate;
4070
4071 for(i = l.firstUnused - 1; i >= end; --i)
4072 {
4073 k = lbeg[i];
4074 val = &lval[k];
4075 idx = &lidx[k];
4076 x1 = 0;
4077 x2 = 0;
4078
4079 for(j = lbeg[i + 1]; j > k; --j)
4080 {
4081 x1 += vec1[*idx] * (*val);
4082 x2 += vec2[*idx++] * (*val++);
4083 }
4084
4085 vec1[lrow[i]] -= x1;
4086
4087 vec2[lrow[i]] -= x2;
4088 }
4089 }
4090
4091 inline int CLUFactorRational::solveUpdateLeft(Rational* vec, int* nonz, int n)
4092 {
4093 int i, j, k, end;
4094 Rational x, y;
4095 Rational* val;
4096 int* lrow, *lidx, *idx;
4097 int* lbeg;
4098
4099 assert(!l.updateType); /* no Forest-Tomlin Updates! */
4100
4101 VectorRational& lval = l.val;
4102 lidx = l.idx;
4103 lrow = l.row;
4104 lbeg = l.start;
4105
4106 end = l.firstUpdate;
4107
4108 for(i = l.firstUnused - 1; i >= end; --i)
4109 {
4110 k = lbeg[i];
4111 assert(k >= 0 && k < l.val.dim());
4112 val = &lval[k];
4113 idx = &lidx[k];
4114 x = 0;
4115
4116 for(j = lbeg[i + 1]; j > k; --j)
4117 {
4118 assert(*idx >= 0 && *idx < thedim);
4119 x += vec[*idx++] * (*val++);
4120 }
4121
4122 k = lrow[i];
4123
4124 y = vec[k];
4125
4126 if(y == 0)
4127 {
4128 y = -x;
4129
4130 if(y != 0)
4131 {
4132 nonz[n++] = k;
4133 vec[k] = y;
4134 }
4135 }
4136 else
4137 {
4138 y -= x;
4139 //vec[k] = ( y != 0 ) ? y : SOPLEX_MARKER;
4140 vec[k] = y;
4141 }
4142 }
4143
4144 return n;
4145 }
4146
4147 inline void CLUFactorRational::solveLeft(Rational* vec, Rational* rhs)
4148 {
4149
4150 if(!l.updateType) /* no Forest-Tomlin Updates */
4151 {
4152 solveUpdateLeft(rhs);
4153 solveUleft(vec, rhs);
4154 solveLleft(vec);
4155 }
4156 else
4157 {
4158 solveUleft(vec, rhs);
4159 solveLleftForest(vec, 0);
4160 solveLleft(vec);
4161 }
4162 }
4163
4164 inline int CLUFactorRational::solveLeftEps(Rational* vec, Rational* rhs, int* nonz)
4165 {
4166
4167 if(!l.updateType) /* no Forest-Tomlin Updates */
4168 {
4169 solveUpdateLeft(rhs);
4170 solveUleft(vec, rhs);
4171 return solveLleftEps(vec, nonz);
4172 }
4173 else
4174 {
4175 solveUleft(vec, rhs);
4176 solveLleftForest(vec, nonz);
4177 return solveLleftEps(vec, nonz);
4178 }
4179 }
4180
4181 inline int CLUFactorRational::solveLeft2(
4182 Rational* vec1,
4183 int* nonz,
4184 Rational* vec2,
4185 Rational* rhs1,
4186 Rational* rhs2)
4187 {
4188
4189 if(!l.updateType) /* no Forest-Tomlin Updates */
4190 {
4191 solveUpdateLeft2(rhs1, rhs2);
4192 solveUleft2(vec1, rhs1, vec2, rhs2);
4193 solveLleft2(vec1, nonz, vec2);
4194 return 0;
4195 }
4196 else
4197 {
4198 solveUleft2(vec1, rhs1, vec2, rhs2);
4199 solveLleft2forest(vec1, nonz, vec2);
4200 solveLleft2(vec1, nonz, vec2);
4201 return 0;
4202 }
4203 }
4204
4205 inline int CLUFactorRational::solveUleft(Rational* vec, int* vecidx,
4206 Rational* rhs, int* rhsidx, int rhsn)
4207 {
4208 Rational x, y;
4209 int i, j, k, n, r, c;
4210 int* rorig, *corig, *cperm;
4211 int* ridx, *rlen, *rbeg, *idx;
4212 Rational* val;
4213
4214 rorig = row.orig;
4215 corig = col.orig;
4216 cperm = col.perm;
4217
4218 /* move rhsidx to a heap
4219 */
4220
4221 for(i = 0; i < rhsn;)
4222 enQueueMinRat(rhsidx, &i, cperm[rhsidx[i]]);
4223
4224 ridx = u.row.idx;
4225
4226 VectorRational& rval = u.row.val;
4227
4228 rlen = u.row.len;
4229
4230 rbeg = u.row.start;
4231
4232 n = 0;
4233
4234 while(rhsn > 0)
4235 {
4236 i = deQueueMinRat(rhsidx, &rhsn);
4237 assert(i >= 0 && i < thedim);
4238 c = corig[i];
4239 assert(c >= 0 && c < thedim);
4240 x = rhs[c];
4241 rhs[c] = 0;
4242
4243 if(x != 0)
4244 {
4245 r = rorig[i];
4246 assert(r >= 0 && r < thedim);
4247 vecidx[n++] = r;
4248 x *= diag[r];
4249 vec[r] = x;
4250 k = rbeg[r];
4251 assert(k >= 0 && k < u.row.val.dim());
4252 idx = &ridx[k];
4253 val = &rval[k];
4254
4255 for(int m = rlen[r]; m; --m)
4256 {
4257 j = *idx++;
4258 assert(j >= 0 && j < thedim);
4259 y = rhs[j];
4260
4261 if(y == 0)
4262 {
4263 y = -x * (*val++);
4264
4265 if(y != 0)
4266 {
4267 rhs[j] = y;
4268 enQueueMinRat(rhsidx, &rhsn, cperm[j]);
4269 }
4270 }
4271 else
4272 {
4273 y -= x * (*val++);
4274 // rhs[j] = ( y != 0 ) ? y : SOPLEX_MARKER;
4275 rhs[j] = y;
4276 }
4277 }
4278 }
4279 }
4280
4281 return n;
4282 }
4283
4284
4285 inline void CLUFactorRational::solveUleftNoNZ(Rational* vec, Rational* rhs, int* rhsidx, int rhsn)
4286 {
4287 Rational x, y;
4288 int i, j, k, r, c;
4289 int* rorig, *corig, *cperm;
4290 int* ridx, *rlen, *rbeg, *idx;
4291 Rational* val;
4292
4293 rorig = row.orig;
4294 corig = col.orig;
4295 cperm = col.perm;
4296
4297 /* move rhsidx to a heap
4298 */
4299
4300 for(i = 0; i < rhsn;)
4301 enQueueMinRat(rhsidx, &i, cperm[rhsidx[i]]);
4302
4303 ridx = u.row.idx;
4304
4305 VectorRational& rval = u.row.val;
4306
4307 rlen = u.row.len;
4308
4309 rbeg = u.row.start;
4310
4311 while(rhsn > 0)
4312 {
4313 i = deQueueMinRat(rhsidx, &rhsn);
4314 assert(i >= 0 && i < thedim);
4315 c = corig[i];
4316 assert(c >= 0 && c < thedim);
4317 x = rhs[c];
4318 rhs[c] = 0;
4319
4320 if(x != 0)
4321 {
4322 r = rorig[i];
4323 assert(r >= 0 && r < thedim);
4324 x *= diag[r];
4325 vec[r] = x;
4326 k = rbeg[r];
4327 assert(k >= 0 && k < u.row.val.dim());
4328 idx = &ridx[k];
4329 val = &rval[k];
4330
4331 for(int m = rlen[r]; m; --m)
4332 {
4333 j = *idx++;
4334 assert(j >= 0 && j < thedim);
4335 y = rhs[j];
4336
4337 if(y == 0)
4338 {
4339 y = -x * (*val++);
4340
4341 if(y != 0)
4342 {
4343 rhs[j] = y;
4344 enQueueMinRat(rhsidx, &rhsn, cperm[j]);
4345 }
4346 }
4347 else
4348 {
4349 y -= x * (*val++);
4350 // rhs[j] = ( y != 0 ) ? y : SOPLEX_MARKER;
4351 rhs[j] = y;
4352 }
4353 }
4354 }
4355 }
4356 }
4357
4358
4359 inline int CLUFactorRational::solveLleftForest(Rational* vec, int* nonz, int n)
4360 {
4361 int i, j, k, end;
4362 Rational x, y;
4363 Rational* val;
4364 int* idx, *lidx, *lrow, *lbeg;
4365
4366 VectorRational& lval = l.val;
4367 lidx = l.idx;
4368 lrow = l.row;
4369 lbeg = l.start;
4370 end = l.firstUpdate;
4371
4372 for(i = l.firstUnused - 1; i >= end; --i)
4373 {
4374 assert(i >= 0 && i < l.val.dim());
4375
4376 if((x = vec[lrow[i]]) != 0)
4377 {
4378 k = lbeg[i];
4379 assert(k >= 0 && k < l.val.dim());
4380 val = &lval[k];
4381 idx = &lidx[k];
4382
4383 for(j = lbeg[i + 1]; j > k; --j)
4384 {
4385 int m = *idx++;
4386 assert(m >= 0 && m < thedim);
4387 y = vec[m];
4388
4389 if(y == 0)
4390 {
4391 y = -x * (*val++);
4392
4393 if(y != 0)
4394 {
4395 vec[m] = y;
4396 nonz[n++] = m;
4397 }
4398 }
4399 else
4400 {
4401 y -= x * (*val++);
4402
4403 }
4404 }
4405 }
4406 }
4407
4408 return n;
4409 }
4410
4411
4412 inline void CLUFactorRational::solveLleftForestNoNZ(Rational* vec)
4413 {
4414 int i, j, k, end;
4415 Rational x;
4416 Rational* val;
4417 int* idx, *lidx, *lrow, *lbeg;
4418
4419 VectorRational& lval = l.val;
4420 lidx = l.idx;
4421 lrow = l.row;
4422 lbeg = l.start;
4423 end = l.firstUpdate;
4424
4425 for(i = l.firstUnused - 1; i >= end; --i)
4426 {
4427 if((x = vec[lrow[i]]) != 0)
4428 {
4429 assert(i >= 0 && i < l.val.dim());
4430 k = lbeg[i];
4431 assert(k >= 0 && k < l.val.dim());
4432 val = &lval[k];
4433 idx = &lidx[k];
4434
4435 for(j = lbeg[i + 1]; j > k; --j)
4436 {
4437 assert(*idx >= 0 && *idx < thedim);
4438 vec[*idx++] -= x * (*val++);
4439 }
4440 }
4441 }
4442 }
4443
4444
4445 inline int CLUFactorRational::solveLleft(Rational* vec, int* nonz, int rn)
4446 {
4447 int i, j, k, n;
4448 int r;
4449 Rational x, y;
4450 Rational* val;
4451 int* ridx, *idx;
4452 int* rbeg;
4453 int* rorig, *rperm;
4454 int* last;
4455
4456 ridx = l.ridx;
4457 VectorRational& rval = l.rval;
4458 rbeg = l.rbeg;
4459 rorig = l.rorig;
4460 rperm = l.rperm;
4461 n = 0;
4462
4463 i = l.firstUpdate - 1;
4464 #ifndef SOPLEX_WITH_L_ROWS
4465 #pragma warn "Not yet implemented, define SOPLEX_WITH_L_ROWS"
4466 VectorRational& lval = l.val;
4467 int* lidx = l.idx;
4468 int* lrow = l.row;
4469 int* lbeg = l.start;
4470
4471 for(; i >= 0; --i)
4472 {
4473 k = lbeg[i];
4474 val = &lval[k];
4475 idx = &lidx[k];
4476 x = 0;
4477
4478 for(j = lbeg[i + 1]; j > k; --j)
4479 x += vec[*idx++] * (*val++);
4480
4481 vec[lrow[i]] -= x;
4482 }
4483
4484 #else
4485
4486 /* move rhsidx to a heap
4487 */
4488 for(i = 0; i < rn;)
4489 enQueueMaxRat(nonz, &i, rperm[nonz[i]]);
4490
4491 last = nonz + thedim;
4492
4493 while(rn > 0)
4494 {
4495 i = deQueueMaxRat(nonz, &rn);
4496 r = rorig[i];
4497 x = vec[r];
4498
4499 if(x != 0)
4500 {
4501 *(--last) = r;
4502 n++;
4503 k = rbeg[r];
4504 j = rbeg[r + 1] - k;
4505 val = &rval[k];
4506 idx = &ridx[k];
4507
4508 while(j-- > 0)
4509 {
4510 assert(l.rperm[*idx] < i);
4511 int m = *idx++;
4512 y = vec[m];
4513
4514 if(y == 0)
4515 {
4516 y = -x * *val++;
4517
4518 if(y != 0)
4519 {
4520 vec[m] = y;
4521 enQueueMaxRat(nonz, &rn, rperm[m]);
4522 }
4523 }
4524 else
4525 {
4526 y -= x * *val++;
4527 // vec[m] = ( y != 0 ) ? y : SOPLEX_MARKER;
4528 vec[m] = y;
4529 }
4530 }
4531 }
4532 else
4533 vec[r] = 0;
4534 }
4535
4536 for(i = 0; i < n; ++i)
4537 *nonz++ = *last++;
4538
4539 #endif
4540
4541 return n;
4542 }
4543
4544
4545 inline void CLUFactorRational::solveLleftNoNZ(Rational* vec)
4546 {
4547 int i, j, k;
4548 int r;
4549 Rational x;
4550 Rational* val;
4551 int* ridx, *idx;
4552 int* rbeg;
4553 int* rorig;
4554
4555 ridx = l.ridx;
4556 VectorRational& rval = l.rval;
4557 rbeg = l.rbeg;
4558 rorig = l.rorig;
4559
4560 #ifndef SOPLEX_WITH_L_ROWS
4561 VectorRational& lval = l.val;
4562 int* lidx = l.idx;
4563 int* lrow = l.row;
4564 int* lbeg = l.start;
4565
4566 i = l.firstUpdate - 1;
4567 assert(i < thedim);
4568
4569 for(; i >= 0; --i)
4570 {
4571 k = lbeg[i];
4572 assert(k >= 0 && k < l.val.dim());
4573 val = &lval[k];
4574 idx = &lidx[k];
4575 x = 0;
4576
4577 for(j = lbeg[i + 1]; j > k; --j)
4578 {
4579 assert(*idx >= 0 && *idx < thedim);
4580 x += vec[*idx++] * (*val++);
4581 }
4582
4583 vec[lrow[i]] -= x;
4584 }
4585
4586 #else
4587
4588 for(i = thedim; i--;)
4589 {
4590 r = rorig[i];
4591 x = vec[r];
4592
4593 if(x != 0)
4594 {
4595 k = rbeg[r];
4596 j = rbeg[r + 1] - k;
4597 val = &rval[k];
4598 idx = &ridx[k];
4599
4600 while(j-- > 0)
4601 {
4602 assert(l.rperm[*idx] < i);
4603 vec[*idx++] -= x * *val++;
4604 }
4605 }
4606 }
4607
4608 #endif
4609 }
4610
4611 inline int CLUFactorRational::vSolveLright(Rational* vec, int* ridx, int rn)
4612 {
4613 int i, j, k, n;
4614 int end;
4615 Rational x;
4616 Rational* val;
4617 int* lrow, *lidx, *idx;
4618 int* lbeg;
4619
4620 VectorRational& lval = l.val;
4621 lidx = l.idx;
4622 lrow = l.row;
4623 lbeg = l.start;
4624
4625 end = l.firstUpdate;
4626
4627 for(i = 0; i < end; ++i)
4628 {
4629 x = vec[lrow[i]];
4630
4631 if(x != 0)
4632 {
4633 k = lbeg[i];
4634 idx = &(lidx[k]);
4635 val = &(lval[k]);
4636
4637 for(j = lbeg[i + 1]; j > k; --j)
4638 {
4639 assert(*idx >= 0 && *idx < thedim);
4640 ridx[rn] = n = *idx++;
4641 rn += (vec[n] == 0) ? 1 : 0;
4642 vec[n] -= x * (*val++);
4643 // vec[n] += ( vec[n] == 0 ) ? SOPLEX_MARKER : 0;
4644 }
4645 }
4646 }
4647
4648 if(l.updateType) /* Forest-Tomlin Updates */
4649 {
4650 end = l.firstUnused;
4651
4652 for(; i < end; ++i)
4653 {
4654 x = 0;
4655 k = lbeg[i];
4656 idx = &(lidx[k]);
4657 val = &(lval[k]);
4658
4659 for(j = lbeg[i + 1]; j > k; --j)
4660 {
4661 assert(*idx >= 0 && *idx < thedim);
4662 x += vec[*idx++] * (*val++);
4663 }
4664
4665 ridx[rn] = j = lrow[i];
4666
4667 rn += (vec[j] == 0) ? 1 : 0;
4668 vec[j] -= x;
4669 // vec[j] += ( vec[j] == 0 ) ? SOPLEX_MARKER : 0;
4670 }
4671 }
4672
4673 return rn;
4674 }
4675
4676 inline void CLUFactorRational::vSolveLright2(Rational* vec, int* ridx, int* rnptr,
4677 Rational* vec2, int* ridx2, int* rn2ptr)
4678 {
4679 int i, j, k, n;
4680 int end;
4681 Rational x, y;
4682 Rational x2, y2;
4683 Rational* val;
4684 int* lrow, *lidx, *idx;
4685 int* lbeg;
4686
4687 int rn = *rnptr;
4688 int rn2 = *rn2ptr;
4689
4690 VectorRational& lval = l.val;
4691 lidx = l.idx;
4692 lrow = l.row;
4693 lbeg = l.start;
4694
4695 end = l.firstUpdate;
4696
4697 for(i = 0; i < end; ++i)
4698 {
4699 j = lrow[i];
4700 x2 = vec2[j];
4701 x = vec[j];
4702
4703 if(x != 0)
4704 {
4705 k = lbeg[i];
4706 idx = &(lidx[k]);
4707 val = &(lval[k]);
4708
4709 if(x2 != 0)
4710 {
4711 for(j = lbeg[i + 1]; j > k; --j)
4712 {
4713 assert(*idx >= 0 && *idx < thedim);
4714 ridx[rn] = ridx2[rn2] = n = *idx++;
4715 y = vec[n];
4716 y2 = vec2[n];
4717 rn += (y == 0) ? 1 : 0;
4718 rn2 += (y2 == 0) ? 1 : 0;
4719 y -= x * (*val);
4720 y2 -= x2 * (*val++);
4721 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 );
4722 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 );
4723 vec[n] = y;
4724 vec2[n] = y2;
4725 }
4726 }
4727 else
4728 {
4729 for(j = lbeg[i + 1]; j > k; --j)
4730 {
4731 assert(*idx >= 0 && *idx < thedim);
4732 ridx[rn] = n = *idx++;
4733 y = vec[n];
4734 rn += (y == 0) ? 1 : 0;
4735 y -= x * (*val++);
4736 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 );
4737 vec[n] = y;
4738 }
4739 }
4740 }
4741 else if(x2 != 0)
4742 {
4743 k = lbeg[i];
4744 idx = &(lidx[k]);
4745 val = &(lval[k]);
4746
4747 for(j = lbeg[i + 1]; j > k; --j)
4748 {
4749 assert(*idx >= 0 && *idx < thedim);
4750 ridx2[rn2] = n = *idx++;
4751 y2 = vec2[n];
4752 rn2 += (y2 == 0) ? 1 : 0;
4753 y2 -= x2 * (*val++);
4754 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 );
4755 vec2[n] = y2;
4756 }
4757 }
4758 }
4759
4760 if(l.updateType) /* Forest-Tomlin Updates */
4761 {
4762 end = l.firstUnused;
4763
4764 for(; i < end; ++i)
4765 {
4766 x = x2 = 0;
4767 k = lbeg[i];
4768 idx = &(lidx[k]);
4769 val = &(lval[k]);
4770
4771 for(j = lbeg[i + 1]; j > k; --j)
4772 {
4773 assert(*idx >= 0 && *idx < thedim);
4774 x += vec[*idx] * (*val);
4775 x2 += vec2[*idx++] * (*val++);
4776 }
4777
4778 ridx[rn] = ridx2[rn2] = j = lrow[i];
4779
4780 rn += (vec[j] == 0) ? 1 : 0;
4781 rn2 += (vec2[j] == 0) ? 1 : 0;
4782 vec[j] -= x;
4783 vec2[j] -= x2;
4784 // vec[j] += ( vec[j] == 0 ) ? SOPLEX_MARKER : 0;
4785 // vec2[j] += ( vec2[j] == 0 ) ? SOPLEX_MARKER : 0;
4786 }
4787 }
4788
4789 *rnptr = rn;
4790
4791 *rn2ptr = rn2;
4792 }
4793
4794 inline void CLUFactorRational::vSolveLright3(Rational* vec, int* ridx, int* rnptr,
4795 Rational* vec2, int* ridx2, int* rn2ptr,
4796 Rational* vec3, int* ridx3, int* rn3ptr)
4797 {
4798 int i, j, k, n;
4799 int end;
4800 Rational x, y;
4801 Rational x2, y2;
4802 Rational x3, y3;
4803 Rational* val;
4804 int* lrow, *lidx, *idx;
4805 int* lbeg;
4806
4807 int rn = *rnptr;
4808 int rn2 = *rn2ptr;
4809 int rn3 = *rn3ptr;
4810
4811 VectorRational& lval = l.val;
4812 lidx = l.idx;
4813 lrow = l.row;
4814 lbeg = l.start;
4815
4816 end = l.firstUpdate;
4817
4818 for(i = 0; i < end; ++i)
4819 {
4820 j = lrow[i];
4821 x3 = vec3[j];
4822 x2 = vec2[j];
4823 x = vec[j];
4824
4825 if(x != 0)
4826 {
4827 k = lbeg[i];
4828 idx = &(lidx[k]);
4829 val = &(lval[k]);
4830
4831 if(x2 != 0)
4832 {
4833 if(x3 != 0)
4834 {
4835 // case 1: all three vectors are nonzero at j
4836 for(j = lbeg[i + 1]; j > k; --j)
4837 {
4838 assert(*idx >= 0 && *idx < thedim);
4839 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++;
4840 y = vec[n];
4841 y2 = vec2[n];
4842 y3 = vec3[n];
4843 rn += (y == 0) ? 1 : 0;
4844 rn2 += (y2 == 0) ? 1 : 0;
4845 rn3 += (y3 == 0) ? 1 : 0;
4846 y -= x * (*val);
4847 y2 -= x2 * (*val);
4848 y3 -= x3 * (*val++);
4849 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 );
4850 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 );
4851 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 );
4852 vec[n] = y;
4853 vec2[n] = y2;
4854 vec3[n] = y3;
4855 }
4856 }
4857 else
4858 {
4859 // case 2: 1 and 2 are nonzero at j
4860 for(j = lbeg[i + 1]; j > k; --j)
4861 {
4862 assert(*idx >= 0 && *idx < thedim);
4863 ridx[rn] = ridx2[rn2] = n = *idx++;
4864 y = vec[n];
4865 y2 = vec2[n];
4866 rn += (y == 0) ? 1 : 0;
4867 rn2 += (y2 == 0) ? 1 : 0;
4868 y -= x * (*val);
4869 y2 -= x2 * (*val++);
4870 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 );
4871 //vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 );
4872 vec[n] = y;
4873 vec2[n] = y2;
4874 }
4875 }
4876 }
4877 else if(x3 != 0)
4878 {
4879 // case 3: 1 and 3 are nonzero at j
4880 for(j = lbeg[i + 1]; j > k; --j)
4881 {
4882 assert(*idx >= 0 && *idx < thedim);
4883 ridx[rn] = ridx3[rn3] = n = *idx++;
4884 y = vec[n];
4885 y3 = vec3[n];
4886 rn += (y == 0) ? 1 : 0;
4887 rn3 += (y3 == 0) ? 1 : 0;
4888 y -= x * (*val);
4889 y3 -= x3 * (*val++);
4890 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 );
4891 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 );
4892 vec[n] = y;
4893 vec3[n] = y3;
4894 }
4895 }
4896 else
4897 {
4898 // case 4: only 1 is nonzero at j
4899 for(j = lbeg[i + 1]; j > k; --j)
4900 {
4901 assert(*idx >= 0 && *idx < thedim);
4902 ridx[rn] = n = *idx++;
4903 y = vec[n];
4904 rn += (y == 0) ? 1 : 0;
4905 y -= x * (*val++);
4906 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 );
4907 vec[n] = y;
4908 }
4909 }
4910 }
4911 else if(x2 != 0)
4912 {
4913 k = lbeg[i];
4914 idx = &(lidx[k]);
4915 val = &(lval[k]);
4916
4917 if(x3 != 0)
4918 {
4919 // case 5: 2 and 3 are nonzero at j
4920 for(j = lbeg[i + 1]; j > k; --j)
4921 {
4922 assert(*idx >= 0 && *idx < thedim);
4923 ridx2[rn2] = ridx3[rn3] = n = *idx++;
4924 y2 = vec2[n];
4925 y3 = vec3[n];
4926 rn2 += (y2 == 0) ? 1 : 0;
4927 rn3 += (y3 == 0) ? 1 : 0;
4928 y2 -= x2 * (*val);
4929 y3 -= x3 * (*val++);
4930 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 );
4931 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 );
4932 vec2[n] = y2;
4933 vec3[n] = y3;
4934 }
4935 }
4936 else
4937 {
4938 // case 6: only 2 is nonzero at j
4939 for(j = lbeg[i + 1]; j > k; --j)
4940 {
4941 assert(*idx >= 0 && *idx < thedim);
4942 ridx2[rn2] = n = *idx++;
4943 y2 = vec2[n];
4944 rn2 += (y2 == 0) ? 1 : 0;
4945 y2 -= x2 * (*val++);
4946 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 );
4947 vec2[n] = y2;
4948 }
4949 }
4950 }
4951 else if(x3 != 0)
4952 {
4953 // case 7: only 3 is nonzero at j
4954 k = lbeg[i];
4955 idx = &(lidx[k]);
4956 val = &(lval[k]);
4957
4958 for(j = lbeg[i + 1]; j > k; --j)
4959 {
4960 assert(*idx >= 0 && *idx < thedim);
4961 ridx3[rn3] = n = *idx++;
4962 y3 = vec3[n];
4963 rn3 += (y3 == 0) ? 1 : 0;
4964 y3 -= x3 * (*val++);
4965 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 );
4966 vec3[n] = y3;
4967 }
4968 }
4969 }
4970
4971 if(l.updateType) /* Forest-Tomlin Updates */
4972 {
4973 end = l.firstUnused;
4974
4975 for(; i < end; ++i)
4976 {
4977 x = x2 = x3 = 0;
4978 k = lbeg[i];
4979 idx = &(lidx[k]);
4980 val = &(lval[k]);
4981
4982 for(j = lbeg[i + 1]; j > k; --j)
4983 {
4984 assert(*idx >= 0 && *idx < thedim);
4985 x += vec[*idx] * (*val);
4986 x2 += vec2[*idx] * (*val);
4987 x3 += vec3[*idx++] * (*val++);
4988 }
4989
4990 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i];
4991
4992 rn += (vec[j] == 0) ? 1 : 0;
4993 rn2 += (vec2[j] == 0) ? 1 : 0;
4994 rn3 += (vec3[j] == 0) ? 1 : 0;
4995 vec[j] -= x;
4996 vec2[j] -= x2;
4997 vec3[j] -= x3;
4998 // vec[j] += ( vec[j] == 0 ) ? SOPLEX_MARKER : 0;
4999 // vec2[j] += ( vec2[j] == 0 ) ? SOPLEX_MARKER : 0;
5000 // vec3[j] += ( vec3[j] == 0 ) ? SOPLEX_MARKER : 0;
5001 }
5002 }
5003
5004 *rnptr = rn;
5005
5006 *rn2ptr = rn2;
5007 *rn3ptr = rn3;
5008 }
5009
5010 inline int CLUFactorRational::vSolveUright(Rational* vec, int* vidx,
5011 Rational* rhs, int* ridx, int rn)
5012 {
5013 int i, j, k, r, c, n;
5014 int* rorig, *corig;
5015 int* rperm;
5016 int* cidx, *clen, *cbeg;
5017 Rational x, y;
5018
5019 int* idx;
5020 Rational* val;
5021
5022 rorig = row.orig;
5023 corig = col.orig;
5024 rperm = row.perm;
5025
5026 cidx = u.col.idx;
5027 VectorRational& cval = u.col.val;
5028 clen = u.col.len;
5029 cbeg = u.col.start;
5030
5031 n = 0;
5032
5033 while(rn > 0)
5034 {
5035 /* Find nonzero with highest permuted row index and setup i and r
5036 */
5037 i = deQueueMaxRat(ridx, &rn);
5038 assert(i >= 0 && i < thedim);
5039 r = rorig[i];
5040 assert(r >= 0 && r < thedim);
5041
5042 x = diag[r] * rhs[r];
5043 rhs[r] = 0;
5044
5045 if(x != 0)
5046 {
5047 c = corig[i];
5048 assert(c >= 0 && c < thedim);
5049 vidx[n++] = c;
5050 vec[c] = x;
5051 val = &cval[cbeg[c]];
5052 idx = &cidx[cbeg[c]];
5053 j = clen[c];
5054
5055 while(j-- > 0)
5056 {
5057 assert(*idx >= 0 && *idx < thedim);
5058 k = *idx++;
5059 assert(k >= 0 && k < thedim);
5060 y = rhs[k];
5061
5062 if(y == 0)
5063 {
5064 y = -x * (*val++);
5065
5066 if(y != 0)
5067 {
5068 rhs[k] = y;
5069 enQueueMaxRat(ridx, &rn, rperm[k]);
5070 }
5071 }
5072 else
5073 {
5074 y -= x * (*val++);
5075 // y += ( y == 0 ) ? SOPLEX_MARKER : 0;
5076 rhs[k] = y;
5077 }
5078 }
5079
5080 if(rn > i * verySparseFactor4rightRat)
5081 {
5082 /* continue with dense case */
5083 for(i = *ridx; i >= 0; --i)
5084 {
5085 r = rorig[i];
5086 assert(r >= 0 && r < thedim);
5087 x = diag[r] * rhs[r];
5088 rhs[r] = 0;
5089
5090 if(x != 0)
5091 {
5092 c = corig[i];
5093 assert(c >= 0 && c < thedim);
5094 vidx[n++] = c;
5095 vec[c] = x;
5096 val = &cval[cbeg[c]];
5097 idx = &cidx[cbeg[c]];
5098 j = clen[c];
5099
5100 while(j-- > 0)
5101 {
5102 assert(*idx >= 0 && *idx < thedim);
5103 rhs[*idx++] -= x * (*val++);
5104 }
5105 }
5106 }
5107
5108 break;
5109 }
5110 }
5111 }
5112
5113 return n;
5114 }
5115
5116
5117 inline void CLUFactorRational::vSolveUrightNoNZ(Rational* vec, Rational* rhs, int* ridx, int rn)
5118 {
5119 int i, j, k, r, c;
5120 int* rorig, *corig;
5121 int* rperm;
5122 int* cidx, *clen, *cbeg;
5123 Rational x, y;
5124
5125 int* idx;
5126 Rational* val;
5127
5128 rorig = row.orig;
5129 corig = col.orig;
5130 rperm = row.perm;
5131
5132 cidx = u.col.idx;
5133 VectorRational& cval = u.col.val;
5134 clen = u.col.len;
5135 cbeg = u.col.start;
5136
5137 while(rn > 0)
5138 {
5139 if(rn > *ridx * verySparseFactor4rightRat)
5140 {
5141 /* continue with dense case */
5142 for(i = *ridx; i >= 0; --i)
5143 {
5144 assert(i >= 0 && i < thedim);
5145 r = rorig[i];
5146 assert(r >= 0 && r < thedim);
5147 x = diag[r] * rhs[r];
5148 rhs[r] = 0;
5149
5150 if(x != 0)
5151 {
5152 c = corig[i];
5153 vec[c] = x;
5154 val = &cval[cbeg[c]];
5155 idx = &cidx[cbeg[c]];
5156 j = clen[c];
5157
5158 while(j-- > 0)
5159 {
5160 assert(*idx >= 0 && *idx < thedim);
5161 rhs[*idx++] -= x * (*val++);
5162 }
5163 }
5164 }
5165
5166 break;
5167 }
5168
5169 /* Find nonzero with highest permuted row index and setup i and r
5170 */
5171 i = deQueueMaxRat(ridx, &rn);
5172
5173 assert(i >= 0 && i < thedim);
5174
5175 r = rorig[i];
5176
5177 assert(r >= 0 && r < thedim);
5178
5179 x = diag[r] * rhs[r];
5180
5181 rhs[r] = 0;
5182
5183 if(x != 0)
5184 {
5185 c = corig[i];
5186 vec[c] = x;
5187 val = &cval[cbeg[c]];
5188 idx = &cidx[cbeg[c]];
5189 j = clen[c];
5190
5191 while(j-- > 0)
5192 {
5193 k = *idx++;
5194 assert(k >= 0 && k < thedim);
5195 y = rhs[k];
5196
5197 if(y == 0)
5198 {
5199 y = -x * (*val++);
5200
5201 if(y != 0)
5202 {
5203 rhs[k] = y;
5204 enQueueMaxRat(ridx, &rn, rperm[k]);
5205 }
5206 }
5207 else
5208 {
5209 y -= x * (*val++);
5210 // y += ( y == 0 ) ? SOPLEX_MARKER : 0;
5211 rhs[k] = y;
5212 }
5213 }
5214 }
5215 }
5216 }
5217
5218
5219 inline int CLUFactorRational::vSolveUright2(Rational* vec, int* vidx, Rational* rhs, int* ridx,
5220 int rn,
5221 Rational* vec2, Rational* rhs2, int* ridx2, int rn2)
5222 {
5223 int i, j, k, r, c, n;
5224 int* rorig, *corig;
5225 int* rperm;
5226 int* cidx, *clen, *cbeg;
5227 Rational x, y;
5228 Rational x2, y2;
5229
5230 int* idx;
5231 Rational* val;
5232
5233 rorig = row.orig;
5234 corig = col.orig;
5235 rperm = row.perm;
5236
5237 cidx = u.col.idx;
5238 VectorRational& cval = u.col.val;
5239 clen = u.col.len;
5240 cbeg = u.col.start;
5241
5242 n = 0;
5243
5244 while(rn + rn2 > 0)
5245 {
5246 /* Find nonzero with highest permuted row index and setup i and r
5247 */
5248 if(rn <= 0)
5249 i = deQueueMaxRat(ridx2, &rn2);
5250 else if(rn2 <= 0)
5251 i = deQueueMaxRat(ridx, &rn);
5252 else if(*ridx2 > *ridx)
5253 i = deQueueMaxRat(ridx2, &rn2);
5254 else if(*ridx2 < *ridx)
5255 i = deQueueMaxRat(ridx, &rn);
5256 else
5257 {
5258 (void) deQueueMaxRat(ridx, &rn);
5259 i = deQueueMaxRat(ridx2, &rn2);
5260 }
5261
5262 assert(i >= 0 && i < thedim);
5263
5264 r = rorig[i];
5265 assert(r >= 0 && r < thedim);
5266
5267 x = diag[r] * rhs[r];
5268 x2 = diag[r] * rhs2[r];
5269 rhs[r] = 0;
5270 rhs2[r] = 0;
5271
5272 if(x != 0)
5273 {
5274 c = corig[i];
5275 vidx[n++] = c;
5276 vec[c] = x;
5277 vec2[c] = x2;
5278 val = &cval[cbeg[c]];
5279 idx = &cidx[cbeg[c]];
5280 j = clen[c];
5281
5282 if(x2 != 0)
5283 {
5284 while(j-- > 0)
5285 {
5286 k = *idx++;
5287 assert(k >= 0 && k < thedim);
5288 y2 = rhs2[k];
5289
5290 if(y2 == 0)
5291 {
5292 y2 = -x2 * (*val);
5293
5294 if(y2 != 0)
5295 {
5296 rhs2[k] = y2;
5297 enQueueMaxRat(ridx2, &rn2, rperm[k]);
5298 }
5299 }
5300 else
5301 {
5302 y2 -= x2 * (*val);
5303 // rhs2[k] = ( y2 != 0 ) ? y2 : SOPLEX_MARKER;
5304 rhs2[k] = y2;
5305 }
5306
5307 y = rhs[k];
5308
5309 if(y == 0)
5310 {
5311 y = -x * (*val++);
5312
5313 if(y != 0)
5314 {
5315 rhs[k] = y;
5316 enQueueMaxRat(ridx, &rn, rperm[k]);
5317 }
5318 }
5319 else
5320 {
5321 y -= x * (*val++);
5322 // y += ( y == 0 ) ? SOPLEX_MARKER : 0;
5323 rhs[k] = y;
5324 }
5325 }
5326 }
5327 else
5328 {
5329 while(j-- > 0)
5330 {
5331 k = *idx++;
5332 assert(k >= 0 && k < thedim);
5333 y = rhs[k];
5334
5335 if(y == 0)
5336 {
5337 y = -x * (*val++);
5338
5339 if(y != 0)
5340 {
5341 rhs[k] = y;
5342 enQueueMaxRat(ridx, &rn, rperm[k]);
5343 }
5344 }
5345 else
5346 {
5347 y -= x * (*val++);
5348 // y += ( y == 0 ) ? SOPLEX_MARKER : 0;
5349 rhs[k] = y;
5350 }
5351 }
5352 }
5353 }
5354 else if(x2 != 0)
5355 {
5356 c = corig[i];
5357 assert(c >= 0 && c < thedim);
5358 vec2[c] = x2;
5359 val = &cval[cbeg[c]];
5360 idx = &cidx[cbeg[c]];
5361 j = clen[c];
5362
5363 while(j-- > 0)
5364 {
5365 k = *idx++;
5366 assert(k >= 0 && k < thedim);
5367 y2 = rhs2[k];
5368
5369 if(y2 == 0)
5370 {
5371 y2 = -x2 * (*val++);
5372
5373 if(y2 != 0)
5374 {
5375 rhs2[k] = y2;
5376 enQueueMaxRat(ridx2, &rn2, rperm[k]);
5377 }
5378 }
5379 else
5380 {
5381 y2 -= x2 * (*val++);
5382 // rhs2[k] = ( y2 != 0 ) ? y2 : SOPLEX_MARKER;
5383 rhs2[k] = y2;
5384 }
5385 }
5386 }
5387
5388 if(rn + rn2 > i * verySparseFactor4rightRat)
5389 {
5390 /* continue with dense case */
5391 if(*ridx > *ridx2)
5392 i = *ridx;
5393 else
5394 i = *ridx2;
5395
5396 for(; i >= 0; --i)
5397 {
5398 assert(i < thedim);
5399 r = rorig[i];
5400 assert(r >= 0 && r < thedim);
5401 x = diag[r] * rhs[r];
5402 x2 = diag[r] * rhs2[r];
5403 rhs[r] = 0;
5404 rhs2[r] = 0;
5405
5406 if(x2 != 0)
5407 {
5408 c = corig[i];
5409 assert(c >= 0 && c < thedim);
5410 vec2[c] = x2;
5411 val = &cval[cbeg[c]];
5412 idx = &cidx[cbeg[c]];
5413 j = clen[c];
5414
5415 if(x != 0)
5416 {
5417 vidx[n++] = c;
5418 vec[c] = x;
5419
5420 while(j-- > 0)
5421 {
5422 assert(*idx >= 0 && *idx < thedim);
5423 rhs[*idx] -= x * (*val);
5424 rhs2[*idx++] -= x2 * (*val++);
5425 }
5426 }
5427 else
5428 {
5429 while(j-- > 0)
5430 {
5431 assert(*idx >= 0 && *idx < thedim);
5432 rhs2[*idx++] -= x2 * (*val++);
5433 }
5434 }
5435 }
5436 else if(x != 0)
5437 {
5438 c = corig[i];
5439 assert(c >= 0 && c < thedim);
5440 vidx[n++] = c;
5441 vec[c] = x;
5442 val = &cval[cbeg[c]];
5443 idx = &cidx[cbeg[c]];
5444 j = clen[c];
5445
5446 while(j-- > 0)
5447 {
5448 assert(*idx >= 0 && *idx < thedim);
5449 rhs[*idx++] -= x * (*val++);
5450 }
5451 }
5452 }
5453
5454 break;
5455 }
5456 }
5457
5458 return n;
5459 }
5460
5461 inline int CLUFactorRational::vSolveUpdateRight(Rational* vec, int* ridx, int n)
5462 {
5463 int i, j, k;
5464 int end;
5465 Rational x, y;
5466 Rational* val;
5467 int* lrow, *lidx, *idx;
5468 int* lbeg;
5469
5470 assert(!l.updateType); /* no Forest-Tomlin Updates */
5471
5472 VectorRational& lval = l.val;
5473 lidx = l.idx;
5474 lrow = l.row;
5475 lbeg = l.start;
5476 end = l.firstUnused;
5477
5478 for(i = l.firstUpdate; i < end; ++i)
5479 {
5480 assert(i >= 0 && i < thedim);
5481 x = vec[lrow[i]];
5482
5483 if(x != 0)
5484 {
5485 k = lbeg[i];
5486 assert(k >= 0 && k < l.val.dim());
5487 idx = &(lidx[k]);
5488 val = &(lval[k]);
5489
5490 for(j = lbeg[i + 1]; j > k; --j)
5491 {
5492 int m = ridx[n] = *idx++;
5493 assert(m >= 0 && m < thedim);
5494 y = vec[m];
5495 n += (y == 0) ? 1 : 0;
5496 y = y - x * (*val++);
5497 // vec[m] = ( y != 0 ) ? y : SOPLEX_MARKER;
5498 vec[m] = y;
5499 }
5500 }
5501 }
5502
5503 return n;
5504 }
5505
5506 inline void CLUFactorRational::vSolveUpdateRightNoNZ(Rational* vec)
5507 {
5508 int i, j, k;
5509 int end;
5510 Rational x;
5511 Rational* val;
5512 int* lrow, *lidx, *idx;
5513 int* lbeg;
5514
5515 assert(!l.updateType); /* no Forest-Tomlin Updates */
5516
5517 VectorRational& lval = l.val;
5518 lidx = l.idx;
5519 lrow = l.row;
5520 lbeg = l.start;
5521 end = l.firstUnused;
5522
5523 for(i = l.firstUpdate; i < end; ++i)
5524 {
5525 assert(i >= 0 && i < thedim);
5526
5527 if((x = vec[lrow[i]]) != 0)
5528 {
5529 k = lbeg[i];
5530 assert(k >= 0 && k < l.val.dim());
5531 idx = &(lidx[k]);
5532 val = &(lval[k]);
5533
5534 for(j = lbeg[i + 1]; j > k; --j)
5535 {
5536 assert(*idx >= 0 && *idx < thedim);
5537 vec[*idx++] -= x * (*val++);
5538 }
5539 }
5540 }
5541 }
5542
5543
5544 inline int CLUFactorRational::vSolveRight4update(Rational* vec,
5545 int* idx, /* result */
5546 Rational* rhs, int* ridx, int rn, /* rhs */
5547 Rational* forest, int* forestNum, int* forestIdx)
5548 {
5549 rn = vSolveLright(rhs, ridx, rn);
5550
5551 /* turn index list into a heap
5552 */
5553
5554 if(forest)
5555 {
5556 Rational x;
5557 int i, j, k;
5558 int* rperm;
5559 int* it = forestIdx;
5560
5561 rperm = row.perm;
5562
5563 for(i = j = 0; i < rn; ++i)
5564 {
5565 k = ridx[i];
5566 assert(k >= 0 && k < thedim);
5567 x = rhs[k];
5568
5569 if(x != 0)
5570 {
5571 enQueueMaxRat(ridx, &j, rperm[*it++ = k]);
5572 forest[k] = x;
5573 }
5574 else
5575 rhs[k] = 0;
5576 }
5577
5578 *forestNum = rn = j;
5579 }
5580 else
5581 {
5582 Rational x;
5583 int i, j, k;
5584 int* rperm;
5585
5586 rperm = row.perm;
5587
5588 for(i = j = 0; i < rn; ++i)
5589 {
5590 k = ridx[i];
5591 assert(k >= 0 && k < thedim);
5592 x = rhs[k];
5593
5594 if(x != 0)
5595 enQueueMaxRat(ridx, &j, rperm[k]);
5596 else
5597 rhs[k] = 0;
5598 }
5599
5600 rn = j;
5601 }
5602
5603 rn = vSolveUright(vec, idx, rhs, ridx, rn);
5604
5605 if(!l.updateType) /* no Forest-Tomlin Updates */
5606 rn = vSolveUpdateRight(vec, idx, rn);
5607
5608 return rn;
5609 }
5610
5611 inline int CLUFactorRational::vSolveRight4update2(Rational* vec,
5612 int* idx, /* result1 */
5613 Rational* rhs, int* ridx, int rn, /* rhs1 */
5614 Rational* vec2, /* result2 */
5615 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */
5616 Rational* forest, int* forestNum, int* forestIdx)
5617 {
5618 vSolveLright2(rhs, ridx, &rn, rhs2, ridx2, &rn2);
5619
5620 /* turn index list into a heap
5621 */
5622
5623 if(forest)
5624 {
5625 Rational x;
5626 int i, j, k;
5627 int* rperm;
5628 int* it = forestIdx;
5629
5630 rperm = row.perm;
5631
5632 for(i = j = 0; i < rn; ++i)
5633 {
5634 k = ridx[i];
5635 assert(k >= 0 && k < thedim);
5636 x = rhs[k];
5637
5638 if(x != 0)
5639 {
5640 enQueueMaxRat(ridx, &j, rperm[*it++ = k]);
5641 forest[k] = x;
5642 }
5643 else
5644 rhs[k] = 0;
5645 }
5646
5647 *forestNum = rn = j;
5648 }
5649 else
5650 {
5651 Rational x;
5652 int i, j, k;
5653 int* rperm;
5654
5655 rperm = row.perm;
5656
5657 for(i = j = 0; i < rn; ++i)
5658 {
5659 k = ridx[i];
5660 assert(k >= 0 && k < thedim);
5661 x = rhs[k];
5662
5663 if(x != 0)
5664 enQueueMaxRat(ridx, &j, rperm[k]);
5665 else
5666 rhs[k] = 0;
5667 }
5668
5669 rn = j;
5670 }
5671
5672 if(rn2 > thedim * verySparseFactor4rightRat)
5673 {
5674 ridx2[0] = thedim - 1;
5675 /* ridx2[1] = thedim - 2; */
5676 }
5677 else
5678 {
5679 Rational x;
5680 /* Rational maxabs; */
5681 int i, j, k;
5682 int* rperm;
5683
5684 /* maxabs = 1; */
5685 rperm = row.perm;
5686
5687 for(i = j = 0; i < rn2; ++i)
5688 {
5689 k = ridx2[i];
5690 assert(k >= 0 && k < thedim);
5691 x = rhs2[k];
5692
5693 if(x == 0)
5694 {
5695 /* maxabs = (maxabs < -x) ? -x : maxabs; */
5696 enQueueMaxRat(ridx2, &j, rperm[k]);
5697 }
5698 else
5699 rhs2[k] = 0;
5700 }
5701
5702 rn2 = j;
5703
5704 }
5705
5706 rn = vSolveUright(vec, idx, rhs, ridx, rn);
5707
5708 vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2);
5709
5710 /*
5711 * rn = vSolveUright2(vec, idx, rhs, ridx, rn, vec2, rhs2, ridx2, rn2 );
5712 */
5713
5714 if(!l.updateType) /* no Forest-Tomlin Updates */
5715 {
5716 rn = vSolveUpdateRight(vec, idx, rn);
5717 vSolveUpdateRightNoNZ(vec2);
5718 }
5719
5720 return rn;
5721 }
5722
5723 inline int CLUFactorRational::vSolveRight4update3(Rational* vec,
5724 int* idx, /* result1 */
5725 Rational* rhs, int* ridx, int rn, /* rhs1 */
5726 Rational* vec2, /* result2 */
5727 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */
5728 Rational* vec3, /* result3 */
5729 Rational* rhs3, int* ridx3, int rn3, /* rhs3 */
5730 Rational* forest, int* forestNum, int* forestIdx)
5731 {
5732
5733 vSolveLright3(rhs, ridx, &rn, rhs2, ridx2, &rn2, rhs3, ridx3, &rn3);
5734 assert(rn >= 0 && rn <= thedim);
5735 assert(rn2 >= 0 && rn2 <= thedim);
5736 assert(rn3 >= 0 && rn3 <= thedim);
5737
5738 /* turn index list into a heap
5739 */
5740
5741 if(forest)
5742 {
5743 Rational x;
5744 int i, j, k;
5745 int* rperm;
5746 int* it = forestIdx;
5747
5748 rperm = row.perm;
5749
5750 for(i = j = 0; i < rn; ++i)
5751 {
5752 k = ridx[i];
5753 assert(k >= 0 && k < thedim);
5754 x = rhs[k];
5755
5756 if(x != 0)
5757 {
5758 enQueueMaxRat(ridx, &j, rperm[*it++ = k]);
5759 forest[k] = x;
5760 }
5761 else
5762 rhs[k] = 0;
5763 }
5764
5765 *forestNum = rn = j;
5766 }
5767 else
5768 {
5769 Rational x;
5770 int i, j, k;
5771 int* rperm;
5772
5773 rperm = row.perm;
5774
5775 for(i = j = 0; i < rn; ++i)
5776 {
5777 k = ridx[i];
5778 assert(k >= 0 && k < thedim);
5779 x = rhs[k];
5780
5781 if(x != 0)
5782 enQueueMaxRat(ridx, &j, rperm[k]);
5783 else
5784 rhs[k] = 0;
5785 }
5786
5787 rn = j;
5788 }
5789
5790 if(rn2 > thedim * verySparseFactor4rightRat)
5791 {
5792 ridx2[0] = thedim - 1;
5793 }
5794 else
5795 {
5796 Rational x;
5797 int i, j, k;
5798 int* rperm;
5799
5800 rperm = row.perm;
5801
5802 for(i = j = 0; i < rn2; ++i)
5803 {
5804 k = ridx2[i];
5805 assert(k >= 0 && k < thedim);
5806 x = rhs2[k];
5807
5808 if(x == 0)
5809 {
5810 enQueueMaxRat(ridx2, &j, rperm[k]);
5811 }
5812 else
5813 rhs2[k] = 0;
5814 }
5815
5816 rn2 = j;
5817 }
5818
5819 if(rn3 > thedim * verySparseFactor4rightRat)
5820 {
5821 ridx3[0] = thedim - 1;
5822 }
5823 else
5824 {
5825 Rational x;
5826 int i, j, k;
5827 int* rperm;
5828
5829 rperm = row.perm;
5830
5831 for(i = j = 0; i < rn3; ++i)
5832 {
5833 k = ridx3[i];
5834 assert(k >= 0 && k < thedim);
5835 x = rhs3[k];
5836
5837 if(x == 0)
5838 {
5839 enQueueMaxRat(ridx3, &j, rperm[k]);
5840 }
5841 else
5842 rhs3[k] = 0;
5843 }
5844
5845 rn3 = j;
5846 }
5847
5848 rn = vSolveUright(vec, idx, rhs, ridx, rn);
5849
5850 vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2);
5851 vSolveUrightNoNZ(vec3, rhs3, ridx3, rn3);
5852
5853 if(!l.updateType) /* no Forest-Tomlin Updates */
5854 {
5855 rn = vSolveUpdateRight(vec, idx, rn);
5856 vSolveUpdateRightNoNZ(vec2);
5857 vSolveUpdateRightNoNZ(vec3);
5858 }
5859
5860 return rn;
5861 }
5862
5863 inline void CLUFactorRational::vSolveRightNoNZ(Rational*
5864 vec2, /* result2 */
5865 Rational* rhs2, int* ridx2, int rn2) /* rhs2 */
5866 {
5867 rn2 = vSolveLright(rhs2, ridx2, rn2);
5868 assert(rn2 >= 0 && rn2 <= thedim);
5869
5870 if(rn2 > thedim * verySparseFactor4rightRat)
5871 {
5872 *ridx2 = thedim - 1;
5873 }
5874 else
5875 {
5876 Rational x;
5877 /* Rational maxabs; */
5878 int i, j, k;
5879 int* rperm;
5880
5881 /* maxabs = 1; */
5882 rperm = row.perm;
5883
5884 for(i = j = 0; i < rn2; ++i)
5885 {
5886 k = ridx2[i];
5887 assert(k >= 0 && k < thedim);
5888 x = rhs2[k];
5889
5890 if(x == 0)
5891 {
5892 /* maxabs = (maxabs < -x) ? -x : maxabs; */
5893 enQueueMaxRat(ridx2, &j, rperm[k]);
5894 }
5895 else
5896 rhs2[k] = 0;
5897 }
5898
5899 rn2 = j;
5900 }
5901
5902 vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2);
5903
5904 if(!l.updateType) /* no Forest-Tomlin Updates */
5905 vSolveUpdateRightNoNZ(vec2);
5906 }
5907
5908 inline int CLUFactorRational::vSolveLeft(Rational* vec, int* idx, /* result */
5909 Rational* rhs, int* ridx, int rn) /* rhs */
5910 {
5911
5912 if(!l.updateType) /* no Forest-Tomlin Updates */
5913 {
5914 rn = solveUpdateLeft(rhs, ridx, rn);
5915 rn = solveUleft(vec, idx, rhs, ridx, rn);
5916 }
5917 else
5918 {
5919 rn = solveUleft(vec, idx, rhs, ridx, rn);
5920 rn = solveLleftForest(vec, idx, rn);
5921 }
5922
5923 return solveLleft(vec, idx, rn);
5924 }
5925
5926 inline int CLUFactorRational::vSolveLeft2(Rational* vec,
5927 int* idx, /* result */
5928 Rational* rhs, int* ridx, int rn, /* rhs */
5929 Rational* vec2, /* result2 */
5930 Rational* rhs2, int* ridx2, int rn2) /* rhs2 */
5931 {
5932
5933 if(!l.updateType) /* no Forest-Tomlin Updates */
5934 {
5935 rn = solveUpdateLeft(rhs, ridx, rn);
5936 rn = solveUleft(vec, idx, rhs, ridx, rn);
5937 rn2 = solveUpdateLeft(rhs2, ridx2, rn2);
5938 solveUleftNoNZ(vec2, rhs2, ridx2, rn2);
5939 }
5940 else
5941 {
5942 rn = solveUleft(vec, idx, rhs, ridx, rn);
5943 rn = solveLleftForest(vec, idx, rn);
5944 solveUleftNoNZ(vec2, rhs2, ridx2, rn2);
5945 solveLleftForestNoNZ(vec2);
5946 }
5947
5948 rn = solveLleft(vec, idx, rn);
5949
5950 solveLleftNoNZ(vec2);
5951
5952 return rn;
5953 }
5954
5955 inline int CLUFactorRational::vSolveLeft3(Rational* vec,
5956 int* idx, /* result */
5957 Rational* rhs, int* ridx, int rn, /* rhs */
5958 Rational* vec2, /* result2 */
5959 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */
5960 Rational* vec3, /* result3 */
5961 Rational* rhs3, int* ridx3, int rn3) /* rhs3 */
5962 {
5963
5964 if(!l.updateType) /* no Forest-Tomlin Updates */
5965 {
5966 rn = solveUpdateLeft(rhs, ridx, rn);
5967 rn = solveUleft(vec, idx, rhs, ridx, rn);
5968 rn2 = solveUpdateLeft(rhs2, ridx2, rn2);
5969 solveUleftNoNZ(vec2, rhs2, ridx2, rn2);
5970 rn3 = solveUpdateLeft(rhs3, ridx3, rn3);
5971 solveUleftNoNZ(vec3, rhs3, ridx3, rn3);
5972 }
5973 else
5974 {
5975 rn = solveUleft(vec, idx, rhs, ridx, rn);
5976 rn = solveLleftForest(vec, idx, rn);
5977 solveUleftNoNZ(vec2, rhs2, ridx2, rn2);
5978 solveLleftForestNoNZ(vec2);
5979 solveUleftNoNZ(vec3, rhs3, ridx3, rn3);
5980 solveLleftForestNoNZ(vec3);
5981 }
5982
5983 rn = solveLleft(vec, idx, rn);
5984
5985 solveLleftNoNZ(vec2);
5986 solveLleftNoNZ(vec3);
5987
5988 return rn;
5989 }
5990
5991 inline void CLUFactorRational::vSolveLeftNoNZ(Rational*
5992 vec2, /* result2 */
5993 Rational* rhs2, int* ridx2, int rn2) /* rhs2 */
5994 {
5995
5996 if(!l.updateType) /* no Forest-Tomlin Updates */
5997 {
5998 rn2 = solveUpdateLeft(rhs2, ridx2, rn2);
5999 solveUleftNoNZ(vec2, rhs2, ridx2, rn2);
6000 }
6001 else
6002 {
6003 solveUleftNoNZ(vec2, rhs2, ridx2, rn2);
6004 solveLleftForestNoNZ(vec2);
6005 }
6006
6007 solveLleftNoNZ(vec2);
6008 }
6009 } // namespace soplex
6010