Scippy

SoPlex

Sequential object-oriented simPlex

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