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