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   	/**@file slufactor.hpp
26   	 * @todo SLUfactor seems to be partly an wrapper for CLUFactor (was C).
27   	 *       This should be properly integrated and demangled.
28   	 * @todo Does is make sense, to call x.clear() when next x.altValues()
29   	 *       is called.
30   	 */
31   	
32   	#include <assert.h>
33   	#include <sstream>
34   	
35   	#ifdef SOPLEX_DEBUG
36   	#include <stdio.h>
37   	#endif
38   	
39   	namespace soplex
40   	{
41   	/// note: we keep this constant since it is just a tradeoff between sparsity and stability and does
42   	/// not need to be changed when precisions are decreased
43   	#define SOPLEX_MINSTABILITY    R(4e-2)
44   	
45   	template <class R>
46   	void SLUFactor<R>::solveRight(VectorBase<R>& x, const VectorBase<R>& b) //const
47   	{
48   	
49   	   this->solveTime->start();
50   	
51   	   this->vec = b;
52   	   x.clear();
53   	   CLUFactor<R>::solveRight(x.get_ptr(), vec.get_ptr());
54   	
55   	   solveCount++;
56   	   solveTime->stop();
57   	}
58   	
59   	template <class R>
60   	void SLUFactor<R>::solveRight(SSVectorBase<R>& x, const SVectorBase<R>& b)  //const
61   	{
62   	
63   	   solveTime->start();
64   	
65   	   vec.assign(b);
66   	   x.clear();
67   	   CLUFactor<R>::solveRight(x.altValues(), vec.get_ptr());
68   	
69   	   solveCount++;
70   	   solveTime->stop();
71   	}
72   	
73   	template <class R>
74   	void SLUFactor<R>::solveRight4update(SSVectorBase<R>& x, const SVectorBase<R>& b)
75   	{
76   	
77   	   solveTime->start();
78   	
79   	   int m;
80   	   int n;
81   	   int f;
82   	
83   	   x.clear();
84   	   ssvec = b;
85   	   n = ssvec.size();
86   	   R epsilon = this->tolerances()->epsilon();
87   	
88   	   if(this->l.updateType == ETA)
89   	   {
90   	      m = this->vSolveRight4update(epsilon, x.altValues(), x.altIndexMem(),
91   	                                   ssvec.altValues(), ssvec.altIndexMem(), n, 0, 0, 0);
92   	      x.setSize(m);
93   	      //x.forceSetup();
94   	      x.unSetup();
95   	      eta.setup_and_assign(x);
96   	   }
97   	   else
98   	   {
99   	      forest.clear();
100  	      m = this->vSolveRight4update(epsilon, x.altValues(), x.altIndexMem(),
101  	                                   ssvec.altValues(), ssvec.altIndexMem(), n,
102  	                                   forest.altValues(), &f, forest.altIndexMem());
103  	      forest.setSize(f);
104  	      forest.forceSetup();
105  	      x.setSize(m);
106  	      x.forceSetup();
107  	   }
108  	
109  	   usetup = true;
110  	   ssvec.setSize(0);
111  	   ssvec.forceSetup();
112  	
113  	   solveCount++;
114  	   solveTime->stop();
115  	}
116  	
117  	template <class R>
118  	void SLUFactor<R>::solve2right4update(
119  	   SSVectorBase<R>&       x,
120  	   VectorBase<R>&        y,
121  	   const SVectorBase<R>& b,
122  	   SSVectorBase<R>&       rhs)
123  	{
124  	
125  	   solveTime->start();
126  	
127  	   int  m;
128  	   int  n;
129  	   int  f;
130  	   int* sidx = ssvec.altIndexMem();
131  	   ssvec.setSize(0);
132  	   ssvec.forceSetup();
133  	   int  rsize = rhs.size();
134  	   int* ridx = rhs.altIndexMem();
135  	   R epsilon = this->tolerances()->epsilon();
136  	
137  	   x.clear();
138  	   y.clear();
139  	   usetup = true;
140  	   ssvec = b;
141  	
142  	   if(this->l.updateType == ETA)
143  	   {
144  	      n = ssvec.size();
145  	      m = this->vSolveRight4update2(epsilon, x.altValues(), x.altIndexMem(),
146  	                                    ssvec.get_ptr(), sidx, n, y.get_ptr(),
147  	                                    epsilon, rhs.altValues(), ridx, rsize, 0, 0, 0);
148  	      x.setSize(m);
149  	      //      x.forceSetup();
150  	      x.unSetup();
151  	      eta.setup_and_assign(x);
152  	   }
153  	   else
154  	   {
155  	      forest.clear();
156  	      n = ssvec.size();
157  	      m = this->vSolveRight4update2(epsilon, x.altValues(), x.altIndexMem(),
158  	                                    ssvec.get_ptr(), sidx, n, y.get_ptr(),
159  	                                    epsilon, rhs.altValues(), ridx, rsize,
160  	                                    forest.altValues(), &f, forest.altIndexMem());
161  	      x.setSize(m);
162  	      x.forceSetup();
163  	      forest.setSize(f);
164  	      forest.forceSetup();
165  	   }
166  	
167  	   rhs.forceSetup();
168  	   ssvec.setSize(0);
169  	   ssvec.forceSetup();
170  	
171  	   solveCount += 2;
172  	   solveTime->stop();
173  	}
174  	
175  	template <class R>
176  	void SLUFactor<R>::solve2right4update(
177  	   SSVectorBase<R>&       x,
178  	   SSVectorBase<R>&       y,
179  	   const SVectorBase<R>& b,
180  	   SSVectorBase<R>&       rhs)
181  	{
182  	
183  	   solveTime->start();
184  	
185  	   int  n;
186  	   int  f;
187  	   int* sidx = ssvec.altIndexMem();
188  	   ssvec.setSize(0);
189  	   ssvec.forceSetup();
190  	   int  rsize = rhs.size();
191  	   int* ridx = rhs.altIndexMem();
192  	   R epsilon = this->tolerances()->epsilon();
193  	
194  	   x.clear();
195  	   y.clear();
196  	   usetup = true;
197  	   ssvec = b;
198  	
199  	   if(this->l.updateType == ETA)
200  	   {
201  	      n = ssvec.size();
202  	      this->vSolveRight4update2sparse(epsilon, x.altValues(), x.altIndexMem(),
203  	                                      ssvec.get_ptr(), sidx, n,
204  	                                      epsilon, y.altValues(), y.altIndexMem(),
205  	                                      rhs.altValues(), ridx, rsize,
206  	                                      0, 0, 0);
207  	      x.setSize(n);
208  	      //      x.forceSetup();
209  	      x.unSetup();
210  	      y.setSize(rsize);
211  	      y.unSetup();
212  	      eta.setup_and_assign(x);
213  	   }
214  	   else
215  	   {
216  	      forest.clear();
217  	      n = ssvec.size();
218  	      this->vSolveRight4update2sparse(epsilon, x.altValues(), x.altIndexMem(),
219  	                                      ssvec.get_ptr(), sidx, n,
220  	                                      epsilon, y.altValues(), y.altIndexMem(),
221  	                                      rhs.altValues(), ridx, rsize,
222  	                                      forest.altValues(), &f, forest.altIndexMem());
223  	      x.setSize(n);
224  	      x.forceSetup();
225  	      y.setSize(rsize);
226  	      y.forceSetup();
227  	      forest.setSize(f);
228  	      forest.forceSetup();
229  	   }
230  	
231  	   rhs.forceSetup();
232  	   ssvec.setSize(0);
233  	   ssvec.forceSetup();
234  	
235  	   solveCount += 2;
236  	   solveTime->stop();
237  	}
238  	
239  	
240  	template <class R>
241  	void SLUFactor<R>::solve3right4update(
242  	   SSVectorBase<R>&       x,
243  	   VectorBase<R>&        y,
244  	   VectorBase<R>&        y2,
245  	   const SVectorBase<R>& b,
246  	   SSVectorBase<R>&       rhs,
247  	   SSVectorBase<R>&       rhs2)
248  	{
249  	
250  	   solveTime->start();
251  	
252  	   int  m;
253  	   int  n;
254  	   int  f;
255  	   int* sidx = ssvec.altIndexMem();
256  	   ssvec.setSize(0);
257  	   ssvec.forceSetup();
258  	   int  rsize = rhs.size();
259  	   int* ridx = rhs.altIndexMem();
260  	   int  rsize2 = rhs2.size();
261  	   int* ridx2 = rhs2.altIndexMem();
262  	   R epsilon = this->tolerances()->epsilon();
263  	
264  	   x.clear();
265  	   y.clear();
266  	   y2.clear();
267  	   usetup = true;
268  	   ssvec = b;
269  	
270  	   if(this->l.updateType == ETA)
271  	   {
272  	      n = ssvec.size();
273  	      m = this->vSolveRight4update3(epsilon,
274  	                                    x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n,
275  	                                    y.get_ptr(), epsilon, rhs.altValues(), ridx, rsize,
276  	                                    y2.get_ptr(), epsilon, rhs2.altValues(), ridx2, rsize2,
277  	                                    0, 0, 0);
278  	      x.setSize(m);
279  	      //      x.forceSetup();
280  	      x.unSetup();
281  	      eta.setup_and_assign(x);
282  	   }
283  	   else
284  	   {
285  	      forest.clear();
286  	      n = ssvec.size();
287  	      m = this->vSolveRight4update3(epsilon,
288  	                                    x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n,
289  	                                    y.get_ptr(), epsilon, rhs.altValues(), ridx, rsize,
290  	                                    y2.get_ptr(), epsilon, rhs2.altValues(), ridx2, rsize2,
291  	                                    forest.altValues(), &f, forest.altIndexMem());
292  	      x.setSize(m);
293  	      x.forceSetup();
294  	      forest.setSize(f);
295  	      forest.forceSetup();
296  	   }
297  	
298  	   rhs.forceSetup();
299  	   rhs2.forceSetup();
300  	   ssvec.setSize(0);
301  	   ssvec.forceSetup();
302  	
303  	   solveCount += 3;
304  	   solveTime->stop();
305  	}
306  	
307  	template <class R>
308  	void SLUFactor<R>::solve3right4update(
309  	   SSVectorBase<R>&       x,
310  	   SSVectorBase<R>&       y,
311  	   SSVectorBase<R>&       y2,
312  	   const SVectorBase<R>& b,
313  	   SSVectorBase<R>&       rhs,
314  	   SSVectorBase<R>&       rhs2)
315  	{
316  	   solveTime->start();
317  	
318  	   int  n;
319  	   int  f;
320  	   int* sidx = ssvec.altIndexMem();
321  	   ssvec.setSize(0);
322  	   ssvec.forceSetup();
323  	   int  rsize = rhs.size();
324  	   int* ridx = rhs.altIndexMem();
325  	   int  rsize2 = rhs2.size();
326  	   int* ridx2 = rhs2.altIndexMem();
327  	   R epsilon = this->tolerances()->epsilon();
328  	
329  	   x.clear();
330  	   y.clear();
331  	   y2.clear();
332  	   usetup = true;
333  	   ssvec = b;
334  	
335  	   if(this->l.updateType == ETA)
336  	   {
337  	      n = ssvec.size();
338  	      this->vSolveRight4update3sparse(epsilon, x.altValues(), x.altIndexMem(),
339  	                                      ssvec.get_ptr(), sidx, n,
340  	                                      epsilon, y.altValues(), y.altIndexMem(),
341  	                                      rhs.altValues(), ridx, rsize,
342  	                                      epsilon, y2.altValues(), y2.altIndexMem(),
343  	                                      rhs2.altValues(), ridx2, rsize2,
344  	                                      0, 0, 0);
345  	      x.setSize(n);
346  	      //      x.forceSetup();
347  	      x.unSetup();
348  	      y.setSize(rsize);
349  	      y.unSetup();
350  	      y2.setSize(rsize2);
351  	      y2.unSetup();
352  	      eta.setup_and_assign(x);
353  	   }
354  	   else
355  	   {
356  	      forest.clear();
357  	      n = ssvec.size();
358  	      this->vSolveRight4update3sparse(epsilon, x.altValues(), x.altIndexMem(),
359  	                                      ssvec.get_ptr(), sidx, n,
360  	                                      epsilon, y.altValues(), y.altIndexMem(),
361  	                                      rhs.altValues(), ridx, rsize,
362  	                                      epsilon, y2.altValues(), y2.altIndexMem(),
363  	                                      rhs2.altValues(), ridx2, rsize2,
364  	                                      forest.altValues(), &f, forest.altIndexMem());
365  	      x.setSize(n);
366  	      x.forceSetup();
367  	      y.setSize(rsize);
368  	      y.forceSetup();
369  	      y2.setSize(rsize2);
370  	      y2.forceSetup();
371  	
372  	      forest.setSize(f);
373  	      forest.forceSetup();
374  	   }
375  	
376  	   rhs.forceSetup();
377  	   rhs2.forceSetup();
378  	   ssvec.setSize(0);
379  	   ssvec.forceSetup();
380  	
381  	   solveCount += 3;
382  	   solveTime->stop();
383  	}
384  	
385  	
386  	template <class R>
387  	void SLUFactor<R>::solveLeft(VectorBase<R>& x, const VectorBase<R>& b) //const
388  	{
389  	
390  	   solveTime->start();
391  	
392  	   vec = b;
393  	   x.clear();
394  	   CLUFactor<R>::solveLeft(x.get_ptr(), vec.get_ptr());
395  	
396  	   solveCount++;
397  	   solveTime->stop();
398  	}
399  	
400  	template <class R>
401  	void SLUFactor<R>::solveLeft(SSVectorBase<R>& x, const SVectorBase<R>& b)  //const
402  	{
403  	   R epsilon = this->tolerances()->epsilon();
404  	
405  	   solveTime->start();
406  	
407  	   // copy to SSVec is done to avoid having to deal with the Nonzero datatype
408  	   // TODO change SVec to standard sparse format
409  	   ssvec.assign(b);
410  	
411  	   x.clear();
412  	   int sz = ssvec.size(); // see .altValues()
413  	   int n = this->vSolveLeft(epsilon, x.altValues(), x.altIndexMem(),
414  	                            ssvec.altValues(), ssvec.altIndexMem(), sz);
415  	
416  	   if(n > 0)
417  	   {
418  	      x.setSize(n);
419  	      x.forceSetup();
420  	   }
421  	   else
422  	      x.unSetup();
423  	
424  	   ssvec.setSize(0);
425  	   ssvec.forceSetup();
426  	
427  	   solveCount++;
428  	   solveTime->stop();
429  	}
430  	
431  	template <class R>
432  	void SLUFactor<R>::solveLeft(
433  	   SSVectorBase<R>&       x,
434  	   VectorBase<R>&        y,
435  	   const SVectorBase<R>& rhs1,
436  	   SSVectorBase<R>&       rhs2) //const
437  	{
438  	
439  	   solveTime->start();
440  	
441  	   int   n;
442  	   R* svec = ssvec.altValues();
443  	   int*  sidx = ssvec.altIndexMem();
444  	   int   rn   = rhs2.size();
445  	   int*  ridx = rhs2.altIndexMem();
446  	   R epsilon = this->tolerances()->epsilon();
447  	
448  	   x.clear();
449  	   y.clear();
450  	   ssvec.assign(rhs1);
451  	   n = ssvec.size(); // see altValues();
452  	   n = this->vSolveLeft2(epsilon, x.altValues(), x.altIndexMem(), svec, sidx, n,
453  	                         y.get_ptr(), rhs2.altValues(), ridx, rn);
454  	
455  	   // this will unsetup x
456  	   x.setSize(n);
457  	
458  	   if(n > 0)
459  	      x.forceSetup();
460  	
461  	   ssvec.setSize(0);
462  	   ssvec.forceSetup();
463  	
464  	   solveCount += 2;
465  	   solveTime->stop();
466  	}
467  	
468  	template <class R>
469  	void SLUFactor<R>::solveLeft(
470  	   SSVectorBase<R>&       x,
471  	   SSVectorBase<R>&       y,
472  	   const SVectorBase<R>& rhs1,
473  	   SSVectorBase<R>&       rhs2) //const
474  	{
475  	
476  	   solveTime->start();
477  	
478  	   int   n1, n2;
479  	   R* svec = ssvec.altValues();
480  	   int*  sidx = ssvec.altIndexMem();
481  	   R epsilon = this->tolerances()->epsilon();
482  	
483  	   x.clear();
484  	   y.clear();
485  	   ssvec.assign(rhs1);
486  	   n1 = ssvec.size(); // see altValues();
487  	   n2 = rhs2.size();
488  	
489  	   if(n2 < 10)
490  	   {
491  	      this->vSolveLeft2sparse(epsilon,
492  	                              x.altValues(), x.altIndexMem(),
493  	                              svec, sidx, n1,
494  	                              y.altValues(), y.altIndexMem(),
495  	                              rhs2.altValues(), rhs2.altIndexMem(), n2);
496  	      y.setSize(n2);
497  	
498  	      if(n2 > 0)
499  	         y.forceSetup();
500  	   }
501  	   else
502  	   {
503  	      n1 = this->vSolveLeft2(epsilon, x.altValues(), x.altIndexMem(), svec, sidx, n1,
504  	                             y.altValues(), rhs2.altValues(), rhs2.altIndexMem(), n2);
505  	      //      y.setup();
506  	   }
507  	
508  	   x.setSize(n1);
509  	
510  	   if(n1 > 0)
511  	      x.forceSetup();
512  	
513  	   //   if (n2 > 0)
514  	   //      y.forceSetup();
515  	
516  	   ssvec.setSize(0);
517  	   ssvec.forceSetup();
518  	
519  	   solveCount += 2;
520  	   solveTime->stop();
521  	}
522  	
523  	
524  	template <class R>
525  	void SLUFactor<R>::solveLeft(
526  	   SSVectorBase<R>&       x,
527  	   VectorBase<R>&        y,
528  	   VectorBase<R>&        z,
529  	   const SVectorBase<R>& rhs1,
530  	   SSVectorBase<R>&       rhs2,
531  	   SSVectorBase<R>&       rhs3)
532  	{
533  	
534  	   solveTime->start();
535  	
536  	   int   n, n2, n3;
537  	   R* svec = ssvec.altValues();
538  	   int*  sidx = ssvec.altIndexMem();
539  	   R epsilon = this->tolerances()->epsilon();
540  	
541  	   x.clear();
542  	   y.clear();
543  	   z.clear();
544  	   ssvec.assign(rhs1);
545  	   n = ssvec.size(); // see altValues();
546  	   n2 = rhs2.size();
547  	   n3 = rhs3.size();
548  	
549  	   n = this->vSolveLeft3(epsilon, x.altValues(), x.altIndexMem(), svec, sidx, n,
550  	                         y.get_ptr(), rhs2.altValues(), rhs2.altIndexMem(), n2,
551  	                         z.get_ptr(), rhs3.altValues(), rhs3.altIndexMem(), n3);
552  	
553  	   x.setSize(n);
554  	
555  	   if(n > 0)
556  	      x.forceSetup();
557  	
558  	   ssvec.setSize(0);
559  	   ssvec.forceSetup();
560  	
561  	   solveCount += 3;
562  	   solveTime->stop();
563  	}
564  	
565  	template <class R>
566  	void SLUFactor<R>::solveLeft(
567  	   SSVectorBase<R>&       x,
568  	   SSVectorBase<R>&       y,
569  	   SSVectorBase<R>&       z,
570  	   const SVectorBase<R>& rhs1,
571  	   SSVectorBase<R>&       rhs2,
572  	   SSVectorBase<R>&       rhs3)
573  	{
574  	
575  	   solveTime->start();
576  	
577  	   int   n1, n2, n3;
578  	   R* svec = ssvec.altValues();
579  	   int*  sidx = ssvec.altIndexMem();
580  	   R epsilon = this->tolerances()->epsilon();
581  	
582  	   x.clear();
583  	   y.clear();
584  	   z.clear();
585  	   ssvec.assign(rhs1);
586  	   n1 = ssvec.size(); // see altValues();
587  	   n2 = rhs2.size();
588  	   n3 = rhs3.size();
589  	   this->vSolveLeft3sparse(epsilon,
590  	                           x.altValues(), x.altIndexMem(),
591  	                           svec, sidx, n1,
592  	                           y.altValues(), y.altIndexMem(),
593  	                           rhs2.altValues(), rhs2.altIndexMem(), n2,
594  	                           z.altValues(), z.altIndexMem(),
595  	                           rhs3.altValues(), rhs3.altIndexMem(), n3);
596  	   x.setSize(n1);
597  	   y.setSize(n2);
598  	   z.setSize(n3);
599  	
600  	   if(n1 > 0)
601  	      x.forceSetup();
602  	
603  	   if(n2 > 0)
604  	      y.forceSetup();
605  	
606  	   if(n3 > 0)
607  	      z.forceSetup();
608  	
609  	   ssvec.setSize(0);
610  	   ssvec.forceSetup();
611  	
612  	   solveCount += 3;
613  	   solveTime->stop();
614  	}
615  	
616  	
617  	template <class R>
618  	R SLUFactor<R>::stability() const
619  	{
620  	   if(status() != this->OK)
621  	      return 0;
622  	
623  	   if(this->maxabs < this->initMaxabs)
624  	      return 1;
625  	
626  	   assert(this->maxabs != 0.0);
627  	   return this->initMaxabs / this->maxabs;
628  	}
629  	
630  	template <class R>
631  	R SLUFactor<R>::matrixMetric(int type) const
632  	{
633  	   R result = 0.0;
634  	
635  	   // catch corner case of empty matrix
636  	   if(dim() == 0)
637  	      return 1.0;
638  	
639  	   switch(type)
640  	   {
641  	   // compute condition estimate by ratio of max/min of elements on the diagonal
642  	   case 0:
643  	   {
644  	      R mindiag = spxAbs(this->diag[0]);
645  	      R maxdiag = spxAbs(this->diag[0]);
646  	
647  	      for(int i = 1; i < dim(); ++i)
648  	      {
649  	         R absdiag = spxAbs(this->diag[i]);
650  	
651  	         if(absdiag < mindiag)
652  	            mindiag = absdiag;
653  	         else if(absdiag > maxdiag)
654  	            maxdiag = absdiag;
655  	      }
656  	
657  	      result = maxdiag / mindiag;
658  	      break;
659  	   }
660  	
661  	   // compute sum of inverses of all elements on the diagonal
662  	   case 1:
663  	      result = 0.0;
664  	
665  	      for(int i = 0; i < dim(); ++i)
666  	         result += 1.0 / this->diag[i];
667  	
668  	      break;
669  	
670  	   // compute determinant (product of all diagonal elements of U)
671  	   case 2:
672  	      result = 1.0;
673  	
674  	      for(int i = 0; i < dim(); ++i)
675  	         result *= this->diag[i];
676  	
677  	      result = 1.0 / result;
678  	      break;
679  	   }
680  	
681  	   return result;
682  	}
683  	
684  	template <class R>
685  	std::string SLUFactor<R>::statistics() const
686  	{
687  	   std::stringstream s;
688  	   s  << "Factorizations     : " << std::setw(10) << getFactorCount() << std::endl
689  	      << "  Time spent       : " << std::setw(10) << std::fixed << std::setprecision(
690  	         2) << getFactorTime() << std::endl
691  	      << "Solves             : " << std::setw(10) << getSolveCount() << std::endl
692  	      << "  Time spent       : " << std::setw(10) << getSolveTime() << std::endl;
693  	
694  	   return s.str();
695  	}
696  	
697  	template <class R>
698  	void SLUFactor<R>::changeEta(int idx, SSVectorBase<R>& et)
699  	{
700  	
701  	   int es = et.size(); // see altValues()
702  	   this->update(idx, et.altValues(), et.altIndexMem(), es);
703  	   et.setSize(0);
704  	   et.forceSetup();
705  	}
706  	
707  	template <class R>
708  	typename SLUFactor<R>::Status SLUFactor<R>::change(
709  	   int             idx,
710  	   const SVectorBase<R>&  subst,
711  	   const SSVectorBase<R>* e)
712  	{
713  	
714  	   // BH 2005-08-23: The boolean usetup indicates that an "update VectorBase<R>"
715  	   // has been set up. I suppose that SSVectorBase<R>  forest is this
716  	   // update VectorBase<R>, which is set up by solveRight4update() and
717  	   // solve2right4update() in order to optimize the basis update.
718  	
719  	   if(usetup)
720  	   {
721  	      if(this->l.updateType == FOREST_TOMLIN)               // FOREST_TOMLIN updates
722  	      {
723  	         // BH 2005-08-19: The size of a SSVectorBase<R>  is the size of the
724  	         // index set, i.e.  the number of nonzeros which is only
725  	         // defined if the SSVectorBase<R>  is set up.  Since
726  	         // SSVectorBase<R> ::altValues() calls unSetup() the size needs to be
727  	         // stored before the following call.
728  	         int fsize = forest.size(); // see altValues()
729  	         this->forestUpdate(idx, forest.altValues(), fsize, forest.altIndexMem());
730  	         forest.setSize(0);
731  	         forest.forceSetup();
732  	      }
733  	      else
734  	      {
735  	         assert(this->l.updateType == ETA);
736  	         changeEta(idx, eta);
737  	      }
738  	   }
739  	   else if(e != 0)                                    // ETA updates
740  	   {
741  	      this->l.updateType = ETA;
742  	      this->updateNoClear(idx, e->values(), e->indexMem(), e->size());
743  	      this->l.updateType = uptype;
744  	   }
745  	   else if(this->l.updateType == FOREST_TOMLIN)             // FOREST_TOMLIN updates
746  	   {
747  	      assert(0);  // probably this part is never called.
748  	      // forestUpdate() with the last parameter set to NULL should fail.
749  	      forest = subst;
750  	      CLUFactor<R>::solveLright(forest.altValues());
751  	      this->forestUpdate(idx, forest.altValues(), 0, nullptr);
752  	      forest.setSize(0);
753  	      forest.forceSetup();
754  	   }
755  	   else                                               // ETA updates
756  	   {
757  	      assert(this->l.updateType == ETA);
758  	      vec = subst;
759  	      eta.clear();
760  	      CLUFactor<R>::solveRight(eta.altValues(), vec.get_ptr());
761  	      changeEta(idx, eta);
762  	   }
763  	
764  	   usetup = false;
765  	
766  	   SPxOut::debug(this, "DSLUFA01\tupdated\t\tstability = {}\n", stability());
767  	
768  	   return status();
769  	}
770  	
771  	template <class R>
772  	void SLUFactor<R>::clear()
773  	{
774  	
775  	   this->rowMemMult    = 5;          /* factor of minimum Memory * #of nonzeros */
776  	   this->colMemMult    = 5;          /* factor of minimum Memory * #of nonzeros */
777  	   this->lMemMult      = 1;          /* factor of minimum Memory * #of nonzeros */
778  	
779  	   this->l.firstUpdate = 0;
780  	   this->l.firstUnused = 0;
781  	   this->thedim        = 0;
782  	
783  	   usetup        = false;
784  	   this->maxabs        = 1;
785  	   this->initMaxabs    = 1;
786  	   lastThreshold = minThreshold;
787  	   minStability  = SOPLEX_MINSTABILITY;
788  	   this->stat          = this->UNLOADED;
789  	
790  	   vec.clear();
791  	   eta.clear();
792  	   ssvec.clear();
793  	   forest.clear();
794  	
795  	   this->u.row.size    = 100;
796  	   this->u.col.size    = 100;
797  	   this->l.size        = 100;
798  	   this->l.startSize   = 100;
799  	
800  	   if(this->l.ridx)
801  	      spx_free(this->l.ridx);
802  	
803  	   if(this->l.rbeg)
804  	      spx_free(this->l.rbeg);
805  	
806  	   if(this->l.rorig)
807  	      spx_free(this->l.rorig);
808  	
809  	   if(this->l.rperm)
810  	      spx_free(this->l.rperm);
811  	
812  	   if(!this->u.row.val.empty())
813  	      this->u.row.val.clear();
814  	
815  	   if(this->u.row.idx)
816  	      spx_free(this->u.row.idx);
817  	
818  	   if(this->u.col.idx)
819  	      spx_free(this->u.col.idx);
820  	
821  	   if(this->l.val.empty())
822  	      this->l.val.clear();
823  	
824  	   if(this->l.idx)
825  	      spx_free(this->l.idx);
826  	
827  	   if(this->l.start)
828  	      spx_free(this->l.start);
829  	
830  	   if(this->l.row)
831  	      spx_free(this->l.row);
832  	
833  	   // G clear() is used in constructor of SLUFactor<R> so we have to
834  	   // G clean up if anything goes wrong here
835  	   try
836  	   {
837  	      this->u.row.val.resize(this->u.row.size);
838  	      spx_alloc(this->u.row.idx, this->u.row.size);
839  	      spx_alloc(this->u.col.idx, this->u.col.size);
840  	
841  	      this->l.val.resize(this->l.size);
842  	      spx_alloc(this->l.idx,   this->l.size);
843  	      spx_alloc(this->l.start, this->l.startSize);
844  	      spx_alloc(this->l.row,   this->l.startSize);
845  	   }
846  	   catch(const SPxMemoryException& x)
847  	   {
848  	      freeAll();
849  	      throw x;
850  	   }
851  	}
852  	
853  	/** assignment used to implement operator=() and copy constructor.
854  	 *  If this is initialised, freeAll() has to be called before.
855  	 *  Class objects from SLUFactor<R> are not copied here.
856  	 */
857  	template <class R>
858  	void SLUFactor<R>::assign(const SLUFactor<R>& old)
859  	{
860  	   this->spxout = old.spxout;
861  	
862  	   solveTime = TimerFactory::createTimer(old.solveTime->type());
863  	   this->factorTime = TimerFactory::createTimer(old.factorTime->type());
864  	
865  	   // slufactor
866  	   uptype        = old.uptype;
867  	   minThreshold  = old.minThreshold;
868  	   minStability  = old.minStability;
869  	   lastThreshold = old.lastThreshold;
870  	
871  	   // clufactor
872  	   this->stat          = old.stat;
873  	   this->thedim        = old.thedim;
874  	   this->nzCnt         = old.nzCnt;
875  	   this->initMaxabs    = old.initMaxabs;
876  	   this->maxabs        = old.maxabs;
877  	   this->rowMemMult    = old.rowMemMult;
878  	   this->colMemMult    = old.colMemMult;
879  	   this->lMemMult      = old.lMemMult;
880  	
881  	   spx_alloc(this->row.perm, this->thedim);
882  	   spx_alloc(this->row.orig, this->thedim);
883  	   spx_alloc(this->col.perm, this->thedim);
884  	   spx_alloc(this->col.orig, this->thedim);
885  	   // spx_alloc(this->diag,     this->thedim);
886  	   this->diag.reserve(this->thedim); // small performance improvement before copying
887  	
888  	   memcpy(this->row.perm, old.row.perm, (unsigned int)this->thedim * sizeof(*this->row.perm));
889  	   memcpy(this->row.orig, old.row.orig, (unsigned int)this->thedim * sizeof(*this->row.orig));
890  	   memcpy(this->col.perm, old.col.perm, (unsigned int)this->thedim * sizeof(*this->col.perm));
891  	   memcpy(this->col.orig, old.col.orig, (unsigned int)this->thedim * sizeof(*this->col.orig));
892  	
893  	   this->diag = old.diag;
894  	
895  	   this->work = vec.get_ptr();
896  	
897  	   /* setup U
898  	    */
899  	   this->u.row.size = old.u.row.size;
900  	   this->u.row.used = old.u.row.used;
901  	
902  	   spx_alloc(this->u.row.elem,  this->thedim);
903  	   this->u.row.val.reserve(this->u.row.size); // small performance improvement
904  	   spx_alloc(this->u.row.idx,   this->u.row.size);
905  	   spx_alloc(this->u.row.start, this->thedim + 1);
906  	   spx_alloc(this->u.row.len,   this->thedim + 1);
907  	   spx_alloc(this->u.row.max,   this->thedim + 1);
908  	
909  	   memcpy(this->u.row.elem,  old.u.row.elem,
910  	          (unsigned int)this->thedim       * sizeof(*this->u.row.elem));
911  	   this->u.row.val = old.u.row.val;
912  	   memcpy(this->u.row.idx,   old.u.row.idx,
913  	          (unsigned int)this->u.row.size   * sizeof(*this->u.row.idx));
914  	   memcpy(this->u.row.start, old.u.row.start,
915  	          (unsigned int)(this->thedim + 1) * sizeof(*this->u.row.start));
916  	   memcpy(this->u.row.len,   old.u.row.len,
917  	          (unsigned int)(this->thedim + 1) * sizeof(*this->u.row.len));
918  	   memcpy(this->u.row.max,   old.u.row.max,
919  	          (unsigned int)(this->thedim + 1) * sizeof(*this->u.row.max));
920  	
921  	   // need to make row list ok ?
922  	   if(this->thedim > 0 && this->stat == this->OK)
923  	   {
924  	      this->u.row.list.idx = old.u.row.list.idx; // .idx neu
925  	
926  	      const typename CLUFactor<R>::Dring* oring = &old.u.row.list;
927  	      typename CLUFactor<R>::Dring*       ring  = &this->u.row.list;
928  	
929  	      while(oring->next != &old.u.row.list)
930  	      {
931  	         ring->next       = &this->u.row.elem[oring->next->idx];
932  	         ring->next->prev = ring;
933  	         oring            = oring->next;
934  	         ring             = ring->next;
935  	      }
936  	
937  	      ring->next       = &this->u.row.list;
938  	      ring->next->prev = ring;
939  	   }
940  	
941  	   this->u.col.size = old.u.col.size;
942  	   this->u.col.used = old.u.col.used;
943  	
944  	   spx_alloc(this->u.col.elem,  this->thedim);
945  	   spx_alloc(this->u.col.idx,   this->u.col.size);
946  	   spx_alloc(this->u.col.start, this->thedim + 1);
947  	   spx_alloc(this->u.col.len,   this->thedim + 1);
948  	   spx_alloc(this->u.col.max,   this->thedim + 1);
949  	
950  	   if(!old.u.col.val.empty())
951  	   {
952  	      this->u.col.val.reserve(this->u.col.size); // small performance improvement before copying
953  	      this->u.col.val = old.u.col.val;
954  	   }
955  	   else
956  	   {
957  	      this->u.col.val.clear();
958  	   }
959  	
960  	   memcpy(this->u.col.elem,  old.u.col.elem,
961  	          (unsigned int)this->thedim       * sizeof(*this->u.col.elem));
962  	   memcpy(this->u.col.idx,   old.u.col.idx,
963  	          (unsigned int)this->u.col.size   * sizeof(*this->u.col.idx));
964  	   memcpy(this->u.col.start, old.u.col.start,
965  	          (unsigned int)(this->thedim + 1) * sizeof(*this->u.col.start));
966  	   memcpy(this->u.col.len,   old.u.col.len,
967  	          (unsigned int)(this->thedim + 1) * sizeof(*this->u.col.len));
968  	   memcpy(this->u.col.max,   old.u.col.max,
969  	          (unsigned int)(this->thedim + 1) * sizeof(*this->u.col.max));
970  	
971  	   // need to make col list ok
972  	   if(this->thedim > 0 && this->stat == this->OK)
973  	   {
974  	      this->u.col.list.idx = old.u.col.list.idx; // .idx neu
975  	
976  	      const typename CLUFactor<R>::Dring* oring = &old.u.col.list;
977  	      typename CLUFactor<R>::Dring*       ring  = &this->u.col.list;
978  	
979  	      while(oring->next != &old.u.col.list)
980  	      {
981  	         ring->next       = &this->u.col.elem[oring->next->idx];
982  	         ring->next->prev = ring;
983  	         oring            = oring->next;
984  	         ring             = ring->next;
985  	      }
986  	
987  	      ring->next       = &this->u.col.list;
988  	      ring->next->prev = ring;
989  	   }
990  	
991  	   /* Setup L
992  	    */
993  	   this->l.size        = old.l.size;
994  	   this->l.startSize   = old.l.startSize;
995  	   this->l.firstUpdate = old.l.firstUpdate;
996  	   this->l.firstUnused = old.l.firstUnused;
997  	   this->l.updateType  = old.l.updateType;
998  	
999  	   this->l.val.reserve(this->l.size); // small performance improvement for copying
1000 	   spx_alloc(this->l.idx,   this->l.size);
1001 	   spx_alloc(this->l.start, this->l.startSize);
1002 	   spx_alloc(this->l.row,   this->l.startSize);
1003 	
1004 	   this->l.val = old.l.val;
1005 	   memcpy(this->l.idx,   old.l.idx, (unsigned int)this->l.size      * sizeof(*this->l.idx));
1006 	   memcpy(this->l.start, old.l.start, (unsigned int)this->l.startSize * sizeof(*this->l.start));
1007 	   memcpy(this->l.row,   old.l.row, (unsigned int)this->l.startSize * sizeof(*this->l.row));
1008 	
1009 	   if(!this->l.rval.empty())
1010 	   {
1011 	      assert(old.l.ridx  != 0);
1012 	      assert(old.l.rbeg  != 0);
1013 	      assert(old.l.rorig != 0);
1014 	      assert(old.l.rperm != 0);
1015 	
1016 	      int memsize = this->l.start[this->l.firstUpdate];
1017 	
1018 	      this->l.rval.reserve(memsize); // small performance improvement for copying
1019 	      spx_alloc(this->l.ridx,  memsize);
1020 	      spx_alloc(this->l.rbeg,  this->thedim + 1);
1021 	      spx_alloc(this->l.rorig, this->thedim);
1022 	      spx_alloc(this->l.rperm, this->thedim);
1023 	
1024 	      this->l.rval = old.l.rval;
1025 	      memcpy(this->l.ridx,  old.l.ridx, (unsigned int)memsize     * sizeof(*this->l.ridx));
1026 	      memcpy(this->l.rbeg,  old.l.rbeg, (unsigned int)(this->thedim + 1) * sizeof(*this->l.rbeg));
1027 	      memcpy(this->l.rorig, old.l.rorig, (unsigned int)this->thedim      * sizeof(*this->l.rorig));
1028 	      memcpy(this->l.rperm, old.l.rperm, (unsigned int)this->thedim      * sizeof(*this->l.rperm));
1029 	   }
1030 	   else
1031 	   {
1032 	      assert(old.l.ridx  == 0);
1033 	      assert(old.l.rbeg  == 0);
1034 	      assert(old.l.rorig == 0);
1035 	      assert(old.l.rperm == 0);
1036 	
1037 	      this->l.ridx  = 0;
1038 	      this->l.rbeg  = 0;
1039 	      this->l.rorig = 0;
1040 	      this->l.rperm = 0;
1041 	   }
1042 	
1043 	   assert(this->row.perm != 0);
1044 	   assert(this->row.orig != 0);
1045 	   assert(this->col.perm != 0);
1046 	   assert(this->col.orig != 0);
1047 	
1048 	   assert(this->u.row.elem  != 0);
1049 	   assert(this->u.row.idx   != 0);
1050 	   assert(this->u.row.start != 0);
1051 	   assert(this->u.row.len   != 0);
1052 	   assert(this->u.row.max   != 0);
1053 	
1054 	   assert(this->u.col.elem  != 0);
1055 	   assert(this->u.col.idx   != 0);
1056 	   assert(this->u.col.start != 0);
1057 	   assert(this->u.col.len   != 0);
1058 	   assert(this->u.col.max   != 0);
1059 	
1060 	   assert(this->l.idx   != 0);
1061 	   assert(this->l.start != 0);
1062 	   assert(this->l.row   != 0);
1063 	
1064 	}
1065 	
1066 	template <class R>
1067 	SLUFactor<R>& SLUFactor<R>::operator=(const SLUFactor<R>& old)
1068 	{
1069 	
1070 	   if(this != &old)
1071 	   {
1072 	      // we don't need to copy them, because they are temporary vectors
1073 	      vec.clear();
1074 	      ssvec.clear();
1075 	
1076 	      eta    = old.eta;
1077 	      forest = old.forest;
1078 	
1079 	      timerType = old.timerType;
1080 	
1081 	      freeAll();
1082 	
1083 	      try
1084 	      {
1085 	         assign(old);
1086 	      }
1087 	      catch(const SPxMemoryException& x)
1088 	      {
1089 	         freeAll();
1090 	         throw x;
1091 	      }
1092 	
1093 	      assert(isConsistent());
1094 	   }
1095 	
1096 	   return *this;
1097 	}
1098 	
1099 	template <class R>
1100 	SLUFactor<R>::SLUFactor()
1101 	   : vec(1)
1102 	   , ssvec(1)
1103 	   , usetup(false)
1104 	   , uptype(FOREST_TOMLIN)
1105 	   , eta(1)
1106 	   , forest(1)
1107 	   , minThreshold(0.01)
1108 	   , timerType(Timer::USER_TIME)
1109 	{
1110 	   this->row.perm    = 0;
1111 	   this->row.orig    = 0;
1112 	   this->col.perm    = 0;
1113 	   this->col.orig    = 0;
1114 	   this->u.row.elem  = 0;
1115 	   this->u.row.idx   = 0;
1116 	   this->u.row.start = 0;
1117 	   this->u.row.len   = 0;
1118 	   this->u.row.max   = 0;
1119 	   this->u.col.elem  = 0;
1120 	   this->u.col.idx   = 0;
1121 	   this->u.col.start = 0;
1122 	   this->u.col.len   = 0;
1123 	   this->u.col.max   = 0;
1124 	   this->l.idx       = 0;
1125 	   this->l.start     = 0;
1126 	   this->l.row       = 0;
1127 	   this->l.ridx      = 0;
1128 	   this->l.rbeg      = 0;
1129 	   this->l.rorig     = 0;
1130 	   this->l.rperm     = 0;
1131 	
1132 	   this->nzCnt  = 0;
1133 	   this->thedim = 0;
1134 	
1135 	   try
1136 	   {
1137 	      solveTime = TimerFactory::createTimer(timerType);
1138 	      this->factorTime = TimerFactory::createTimer(timerType);
1139 	      spx_alloc(this->row.perm, this->thedim);
1140 	      spx_alloc(this->row.orig, this->thedim);
1141 	      spx_alloc(this->col.perm, this->thedim);
1142 	      spx_alloc(this->col.orig, this->thedim);
1143 	
1144 	      this->diag.resize(this->thedim);
1145 	
1146 	      this->work = vec.get_ptr();
1147 	
1148 	      this->u.row.size = 1;
1149 	      this->u.row.used = 0;
1150 	      spx_alloc(this->u.row.elem,  this->thedim);
1151 	      this->u.row.val.resize(this->u.row.size);
1152 	      spx_alloc(this->u.row.idx,   this->u.row.size);
1153 	      spx_alloc(this->u.row.start, this->thedim + 1);
1154 	      spx_alloc(this->u.row.len,   this->thedim + 1);
1155 	      spx_alloc(this->u.row.max,   this->thedim + 1);
1156 	
1157 	      this->u.row.list.idx      = this->thedim;
1158 	      this->u.row.start[this->thedim] = 0;
1159 	      this->u.row.max  [this->thedim] = 0;
1160 	      this->u.row.len  [this->thedim] = 0;
1161 	
1162 	      this->u.col.size = 1;
1163 	      this->u.col.used = 0;
1164 	      spx_alloc(this->u.col.elem,  this->thedim);
1165 	      spx_alloc(this->u.col.idx,   this->u.col.size);
1166 	      spx_alloc(this->u.col.start, this->thedim + 1);
1167 	      spx_alloc(this->u.col.len,   this->thedim + 1);
1168 	      spx_alloc(this->u.col.max,   this->thedim + 1);
1169 	
1170 	      this->u.col.list.idx      = this->thedim;
1171 	      this->u.col.start[this->thedim] = 0;
1172 	      this->u.col.max[this->thedim]   = 0;
1173 	      this->u.col.len[this->thedim]   = 0;
1174 	
1175 	      this->l.size = 1;
1176 	
1177 	      this->l.val.resize(this->l.size);
1178 	      spx_alloc(this->l.idx, this->l.size);
1179 	
1180 	      this->l.startSize   = 1;
1181 	      this->l.firstUpdate = 0;
1182 	      this->l.firstUnused = 0;
1183 	
1184 	      spx_alloc(this->l.start, this->l.startSize);
1185 	      spx_alloc(this->l.row,   this->l.startSize);
1186 	   }
1187 	   catch(const SPxMemoryException& x)
1188 	   {
1189 	      freeAll();
1190 	      throw x;
1191 	   }
1192 	
1193 	   this->l.ridx  = 0;
1194 	   this->l.rbeg  = 0;
1195 	   this->l.rorig = 0;
1196 	   this->l.rperm = 0;
1197 	
1198 	   SLUFactor<R>::clear(); // clear() is virtual
1199 	
1200 	   this->factorCount = 0;
1201 	   this->hugeValues = 0;
1202 	   solveCount  = 0;
1203 	   assert(this->row.perm != 0);
1204 	
1205 	   assert(this->row.orig != 0);
1206 	   assert(this->col.perm != 0);
1207 	   assert(this->col.orig != 0);
1208 	
1209 	   assert(this->u.row.elem  != 0);
1210 	   assert(this->u.row.idx   != 0);
1211 	   assert(this->u.row.start != 0);
1212 	   assert(this->u.row.len   != 0);
1213 	   assert(this->u.row.max   != 0);
1214 	
1215 	   assert(this->u.col.elem  != 0);
1216 	   assert(this->u.col.idx   != 0);
1217 	   assert(this->u.col.start != 0);
1218 	   assert(this->u.col.len   != 0);
1219 	   assert(this->u.col.max   != 0);
1220 	
1221 	   assert(this->l.idx   != 0);
1222 	   assert(this->l.start != 0);
1223 	   assert(this->l.row   != 0);
1224 	
1225 	   assert(SLUFactor<R>::isConsistent());
1226 	}
1227 	
1228 	template <class R>
1229 	SLUFactor<R>::SLUFactor(const SLUFactor<R>& old)
1230 	   : SLinSolver<R>(old)
1231 	   , vec(1)     // we don't need to copy it, because they are temporary vectors
1232 	   , ssvec(1)   // we don't need to copy it, because they are temporary vectors
1233 	   , usetup(old.usetup)
1234 	   , eta(old.eta)
1235 	   , forest(old.forest)
1236 	   , timerType(old.timerType)
1237 	{
1238 	   this->row.perm    = 0;
1239 	   this->row.orig    = 0;
1240 	   this->col.perm    = 0;
1241 	   this->col.orig    = 0;
1242 	   this->u.row.elem  = 0;
1243 	   this->u.row.val.clear();
1244 	   this->u.row.idx   = 0;
1245 	   this->u.row.start = 0;
1246 	   this->u.row.len   = 0;
1247 	   this->u.row.max   = 0;
1248 	   this->u.col.elem  = 0;
1249 	   this->u.col.idx   = 0;
1250 	   this->u.col.start = 0;
1251 	   this->u.col.len   = 0;
1252 	   this->u.col.max   = 0;
1253 	   this->l.idx       = 0;
1254 	   this->l.start     = 0;
1255 	   this->l.row       = 0;
1256 	   this->l.ridx      = 0;
1257 	   this->l.rbeg      = 0;
1258 	   this->l.rorig     = 0;
1259 	   this->l.rperm     = 0;
1260 	
1261 	   solveCount = 0;
1262 	
1263 	   try
1264 	   {
1265 	      assign(old);
1266 	   }
1267 	   catch(const SPxMemoryException& x)
1268 	   {
1269 	      freeAll();
1270 	      throw x;
1271 	   }
1272 	
1273 	   assert(SLUFactor<R>::isConsistent());
1274 	}
1275 	
1276 	template <class R>
1277 	void SLUFactor<R>::freeAll()
1278 	{
1279 	
1280 	   if(this->row.perm)
1281 	      spx_free(this->row.perm);
1282 	
1283 	   if(this->row.orig)
1284 	      spx_free(this->row.orig);
1285 	
1286 	   if(this->col.perm)
1287 	      spx_free(this->col.perm);
1288 	
1289 	   if(this->col.orig)
1290 	      spx_free(this->col.orig);
1291 	
1292 	   if(this->u.row.elem)
1293 	      spx_free(this->u.row.elem);
1294 	
1295 	   if(!this->u.row.val.empty())
1296 	      this->u.row.val.clear();
1297 	
1298 	   if(this->u.row.idx)
1299 	      spx_free(this->u.row.idx);
1300 	
1301 	   if(this->u.row.start)
1302 	      spx_free(this->u.row.start);
1303 	
1304 	   if(this->u.row.len)
1305 	      spx_free(this->u.row.len);
1306 	
1307 	   if(this->u.row.max)
1308 	      spx_free(this->u.row.max);
1309 	
1310 	   if(this->u.col.elem)
1311 	      spx_free(this->u.col.elem);
1312 	
1313 	   if(this->u.col.idx)
1314 	      spx_free(this->u.col.idx);
1315 	
1316 	   if(this->u.col.start)
1317 	      spx_free(this->u.col.start);
1318 	
1319 	   if(this->u.col.len)
1320 	      spx_free(this->u.col.len);
1321 	
1322 	   if(this->u.col.max)
1323 	      spx_free(this->u.col.max);
1324 	
1325 	   if(!this->l.val.empty())
1326 	      this->l.val.clear();
1327 	
1328 	   if(this->l.idx)
1329 	      spx_free(this->l.idx);
1330 	
1331 	   if(this->l.start)
1332 	      spx_free(this->l.start);
1333 	
1334 	   if(this->l.row)
1335 	      spx_free(this->l.row);
1336 	
1337 	
1338 	   if(!this->u.col.val.empty())
1339 	      this->u.col.val.clear();
1340 	
1341 	   if(this->l.ridx)
1342 	      spx_free(this->l.ridx);
1343 	
1344 	   if(this->l.rbeg)
1345 	      spx_free(this->l.rbeg);
1346 	
1347 	   if(this->l.rorig)
1348 	      spx_free(this->l.rorig);
1349 	
1350 	   if(this->l.rperm)
1351 	      spx_free(this->l.rperm);
1352 	
1353 	   if(solveTime)
1354 	   {
1355 	      solveTime->~Timer();
1356 	      spx_free(solveTime);
1357 	   }
1358 	
1359 	   if(this->factorTime)
1360 	   {
1361 	      this->factorTime->~Timer();
1362 	      spx_free(this->factorTime);
1363 	   }
1364 	}
1365 	
1366 	template <class R>
1367 	SLUFactor<R>::~SLUFactor()
1368 	{
1369 	   freeAll();
1370 	}
1371 	
1372 	template <class R>
1373 	static R betterThreshold(R th, Real epsilon)
1374 	{
1375 	   assert(th < R(1.0));
1376 	
1377 	   if(LT(th, R(0.1), epsilon))
1378 	      th *= R(10.0);
1379 	   else if(LT(th, R(0.9), epsilon))
1380 	      th = (th + R(1.0)) / R(2.0);
1381 	   else if(LT(th, R(0.999), epsilon))
1382 	      th = R(0.99999);
1383 	
1384 	   assert(th < R(1.0));
1385 	
1386 	   return th;
1387 	}
1388 	
1389 	template <class R>
1390 	typename SLUFactor<R>::Status SLUFactor<R>::load(const SVectorBase<R>* matrix[], int dm)
1391 	{
1392 	   assert(dm     >= 0);
1393 	   assert(matrix != 0);
1394 	
1395 	   R lastStability = stability();
1396 	
1397 	   initDR(this->u.row.list);
1398 	   initDR(this->u.col.list);
1399 	
1400 	   usetup        = false;
1401 	   this->l.updateType  = uptype;
1402 	   this->l.firstUpdate = 0;
1403 	   this->l.firstUnused = 0;
1404 	
1405 	   if(dm != this->thedim)
1406 	   {
1407 	      clear();
1408 	
1409 	      this->thedim = dm;
1410 	      vec.reDim(this->thedim);
1411 	      ssvec.reDim(this->thedim);
1412 	      eta.reDim(this->thedim);
1413 	      forest.reDim(this->thedim);
1414 	      this->work = vec.get_ptr();
1415 	
1416 	      spx_realloc(this->row.perm, this->thedim);
1417 	      spx_realloc(this->row.orig, this->thedim);
1418 	      spx_realloc(this->col.perm, this->thedim);
1419 	      spx_realloc(this->col.orig, this->thedim);
1420 	      this->diag.resize(this->thedim);
1421 	
1422 	      spx_realloc(this->u.row.elem,  this->thedim);
1423 	      spx_realloc(this->u.row.len,   this->thedim + 1);
1424 	      spx_realloc(this->u.row.max,   this->thedim + 1);
1425 	      spx_realloc(this->u.row.start, this->thedim + 1);
1426 	
1427 	      spx_realloc(this->u.col.elem,  this->thedim);
1428 	      spx_realloc(this->u.col.len,   this->thedim + 1);
1429 	      spx_realloc(this->u.col.max,   this->thedim + 1);
1430 	      spx_realloc(this->u.col.start, this->thedim + 1);
1431 	
1432 	      this->l.startSize = this->thedim + SOPLEX_MAXUPDATES;
1433 	
1434 	      spx_realloc(this->l.row,   this->l.startSize);
1435 	      spx_realloc(this->l.start, this->l.startSize);
1436 	   }
1437 	   // the last factorization was reasonably stable, so we decrease the Markowitz threshold (stored in lastThreshold) in
1438 	   // order to favour sparsity
1439 	   else if(lastStability > 2.0 * minStability)
1440 	   {
1441 	      // we reset lastThreshold to its previous value in the sequence minThreshold, betterThreshold(minThreshold),
1442 	      // betterThreshold(betterThreshold(minThreshold)), ...
1443 	      R last   = minThreshold;
1444 	      R better = betterThreshold(last, this->tolerances()->epsilon());
1445 	
1446 	      while(better < lastThreshold)
1447 	      {
1448 	         last   = better;
1449 	         better = betterThreshold(last, this->tolerances()->epsilon());
1450 	      }
1451 	
1452 	      lastThreshold = last;
1453 	
1454 	      // we reset the minimum stability (which might have been decreased below) to ensure that the increased sparsity
1455 	      // does not hurt the stability
1456 	      minStability  = 2 * SOPLEX_MINSTABILITY;
1457 	   }
1458 	
1459 	   this->u.row.list.idx      = this->thedim;
1460 	   this->u.row.start[this->thedim] = 0;
1461 	   this->u.row.max[this->thedim]   = 0;
1462 	   this->u.row.len[this->thedim]   = 0;
1463 	
1464 	   this->u.col.list.idx      = this->thedim;
1465 	   this->u.col.start[this->thedim] = 0;
1466 	   this->u.col.max[this->thedim]   = 0;
1467 	   this->u.col.len[this->thedim]   = 0;
1468 	
1469 	   for(;;)
1470 	   {
1471 	      ///@todo if the factorization fails with stat = SINGULAR, distinuish between proven singularity (e.g., because of
1472 	      ///an empty column) and singularity due to numerics, that could be avoided by changing minStability and
1473 	      ///lastThreshold; in the first case, we want to abort, otherwise change the numerics
1474 	      this->stat = this->OK;
1475 	      this->factor(matrix, lastThreshold, this->tolerances()->epsilonFactorization());
1476 	
1477 	      // finish if the factorization is stable
1478 	      if(stability() >= minStability)
1479 	         break;
1480 	
1481 	      // otherwise, we increase the Markowitz threshold
1482 	      R x = lastThreshold;
1483 	      lastThreshold = betterThreshold(lastThreshold, this->tolerances()->epsilon());
1484 	
1485 	      // until it doesn't change anymore at its maximum value
1486 	      if(EQ(x, lastThreshold, this->tolerances()->epsilon()))
1487 	         break;
1488 	
1489 	      // we relax the stability requirement
1490 	      minStability /= 2.0;
1491 	
1492 	      SPX_MSG_INFO3((*this->spxout), (*this->spxout) <<
1493 	                    "ISLUFA01 refactorizing with increased Markowitz threshold: "
1494 	                    << lastThreshold << std::endl;)
1495 	   }
1496 	
1497 	   SPxOut::debug(this, "DSLUFA02 threshold = {} \tstability = {}\tminStability = {}\n", lastThreshold,
1498 	                 stability(), minStability);
1499 	   SPX_DEBUG(
1500 	   {
1501 	      int i;
1502 	      FILE* fl = fopen("dump.lp", "w");
1503 	      std::cout << "DSLUFA03 Basis:\n";
1504 	      int j = 0;
1505 	
1506 	      for(i = 0; i < dim(); ++i)
1507 	         j += matrix[i]->size();
1508 	      for(i = 0; i < dim(); ++i)
1509 	      {
1510 	         for(j = 0; j < matrix[i]->size(); ++j)
1511 	         {
1512 	            fprintf(fl, "%8d  %8d  ",
1513 	                    i + 1, matrix[i]->index(j) + 1);
1514 	            std::cout << matrix[i]->value(j) << std::endl;
1515 	         }
1516 	      }
1517 	      fclose(fl);
1518 	      std::cout << "DSLUFA04 LU-Factors:" << std::endl;
1519 	      dump();
1520 	
1521 	      std::cout << "DSLUFA05 threshold = " << lastThreshold
1522 	                << "\tstability = " << stability() << std::endl;
1523 	   }
1524 	   )
1525 	
1526 	   assert(isConsistent());
1527 	   return Status(this->stat);
1528 	}
1529 	
1530 	
1531 	template <class R>
1532 	bool SLUFactor<R>::isConsistent() const
1533 	{
1534 	#ifdef ENABLE_CONSISTENCY_CHECKS
1535 	   return CLUFactor<R>::isConsistent();
1536 	#else
1537 	   return true;
1538 	#endif
1539 	}
1540 	
1541 	template <class R>
1542 	void SLUFactor<R>::dump() const
1543 	{
1544 	   CLUFactor<R>::dump();
1545 	}
1546 	} // namespace soplex
1547