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