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_rational.hpp
26   	 * @todo SLUFactorRational seems to be partly an wrapper for CLUFactorRational (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   	#include "soplex/cring.h"
35   	#include "soplex/spxdefines.h"
36   	
37   	#ifdef SOPLEX_DEBUG
38   	#include <stdio.h>
39   	#endif
40   	
41   	namespace soplex
42   	{
43   	#define SOPLEX_MINSTABILITYRAT    SOPLEX_REAL(4e-2)
44   	
45   	inline void SLUFactorRational::solveRight(VectorRational& x, const VectorRational& b) //const
46   	{
47   	
48   	   solveTime->start();
49   	
50   	   vec = b;
51   	   CLUFactorRational::solveRight(x.get_ptr(), vec.get_ptr());
52   	
53   	   solveCount++;
54   	   solveTime->stop();
55   	}
56   	
57   	inline void SLUFactorRational::solveRight(SSVectorRational& x, const SVectorRational& b) //const
58   	{
59   	
60   	   solveTime->start();
61   	
62   	   vec.assign(b);
63   	   x.clear();
64   	   CLUFactorRational::solveRight(x.altValues(), vec.get_ptr());
65   	
66   	   solveCount++;
67   	   solveTime->stop();
68   	}
69   	
70   	inline void SLUFactorRational::solveRight4update(SSVectorRational& x, const SVectorRational& b)
71   	{
72   	
73   	   solveTime->start();
74   	
75   	   int m;
76   	   int n;
77   	   int f = 0;
78   	
79   	   x.clear();
80   	   ssvec = b;
81   	   n = ssvec.size();
82   	
83   	   if(l.updateType == ETA)
84   	   {
85   	      m = vSolveRight4update(x.altValues(), x.altIndexMem(),
86   	                             ssvec.altValues(), ssvec.altIndexMem(), n, 0, 0, 0);
87   	      x.setSize(m);
88   	      //x.forceSetup();
89   	      x.unSetup();
90   	      eta.setup_and_assign(x);
91   	   }
92   	   else
93   	   {
94   	      forest.clear();
95   	      m = vSolveRight4update(x.altValues(), x.altIndexMem(),
96   	                             ssvec.altValues(), ssvec.altIndexMem(), n,
97   	                             forest.altValues(), &f, forest.altIndexMem());
98   	      forest.setSize(f);
99   	      forest.forceSetup();
100  	      x.setSize(m);
101  	      x.forceSetup();
102  	   }
103  	
104  	   usetup = true;
105  	
106  	   solveCount++;
107  	   solveTime->stop();
108  	}
109  	
110  	inline void SLUFactorRational::solve2right4update(
111  	   SSVectorRational&      x,
112  	   VectorRational&        y,
113  	   const SVectorRational& b,
114  	   SSVectorRational&      rhs)
115  	{
116  	
117  	   solveTime->start();
118  	
119  	   int  m;
120  	   int  n;
121  	   int  f = 0;
122  	   int* sidx = ssvec.altIndexMem();
123  	   int  rsize = rhs.size();
124  	   int* ridx = rhs.altIndexMem();
125  	
126  	   x.clear();
127  	   y.clear();
128  	   usetup = true;
129  	   ssvec = b;
130  	
131  	   if(l.updateType == ETA)
132  	   {
133  	      n = ssvec.size();
134  	      m = vSolveRight4update2(x.altValues(), x.altIndexMem(),
135  	                              ssvec.get_ptr(), sidx, n, y.get_ptr(),
136  	                              rhs.altValues(), ridx, rsize, 0, 0, 0);
137  	      x.setSize(m);
138  	      //      x.forceSetup();
139  	      x.unSetup();
140  	      eta.setup_and_assign(x);
141  	   }
142  	   else
143  	   {
144  	      forest.clear();
145  	      n = ssvec.size();
146  	      m = vSolveRight4update2(x.altValues(), x.altIndexMem(),
147  	                              ssvec.get_ptr(), sidx, n, y.get_ptr(),
148  	                              rhs.altValues(), ridx, rsize,
149  	                              forest.altValues(), &f, forest.altIndexMem());
150  	      x.setSize(m);
151  	      x.forceSetup();
152  	      forest.setSize(f);
153  	      forest.forceSetup();
154  	   }
155  	
156  	   solveCount++;
157  	   solveTime->stop();
158  	}
159  	
160  	inline void SLUFactorRational::solve3right4update(
161  	   SSVectorRational&      x,
162  	   VectorRational&        y,
163  	   VectorRational&        y2,
164  	   const SVectorRational& b,
165  	   SSVectorRational&      rhs,
166  	   SSVectorRational&      rhs2)
167  	{
168  	
169  	   solveTime->start();
170  	
171  	   int  m;
172  	   int  n;
173  	   int  f;
174  	   int* sidx = ssvec.altIndexMem();
175  	   int  rsize = rhs.size();
176  	   int* ridx = rhs.altIndexMem();
177  	   int rsize2 = rhs2.size();
178  	   int* ridx2 = rhs2.altIndexMem();
179  	
180  	   x.clear();
181  	   y.clear();
182  	   y2.clear();
183  	   usetup = true;
184  	   ssvec = b;
185  	
186  	   if(l.updateType == ETA)
187  	   {
188  	      n = ssvec.size();
189  	      m = vSolveRight4update3(x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n,
190  	                              y.get_ptr(), rhs.altValues(), ridx, rsize,
191  	                              y2.get_ptr(), rhs2.altValues(), ridx2, rsize2,
192  	                              0, 0, 0);
193  	      x.setSize(m);
194  	      //      x.forceSetup();
195  	      x.unSetup();
196  	      eta.setup_and_assign(x);
197  	   }
198  	   else
199  	   {
200  	      forest.clear();
201  	      n = ssvec.size();
202  	      m = vSolveRight4update3(x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n,
203  	                              y.get_ptr(), rhs.altValues(), ridx, rsize,
204  	                              y2.get_ptr(), rhs2.altValues(), ridx2, rsize2,
205  	                              forest.altValues(), &f, forest.altIndexMem());
206  	      x.setSize(m);
207  	      x.forceSetup();
208  	      forest.setSize(f);
209  	      forest.forceSetup();
210  	   }
211  	
212  	   solveCount++;
213  	   solveTime->stop();
214  	}
215  	
216  	inline void SLUFactorRational::solveLeft(VectorRational& x, const VectorRational& b) //const
217  	{
218  	
219  	   solveTime->start();
220  	
221  	   vec = b;
222  	   ///@todo Why is x.clear() here used and not with solveRight() ?
223  	   x.clear();
224  	   CLUFactorRational::solveLeft(x.get_ptr(), vec.get_ptr());
225  	
226  	   solveCount++;
227  	   solveTime->stop();
228  	}
229  	
230  	inline void SLUFactorRational::solveLeft(SSVectorRational& x, const SVectorRational& b) //const
231  	{
232  	
233  	   solveTime->start();
234  	
235  	   ssvec.assign(b);
236  	
237  	   x.clear();
238  	   int sz = ssvec.size(); // see .altValues()
239  	   int n = vSolveLeft(x.altValues(), x.altIndexMem(),
240  	                      ssvec.altValues(), ssvec.altIndexMem(), sz);
241  	
242  	   if(n > 0)
243  	   {
244  	      x.setSize(n);
245  	      x.forceSetup();
246  	   }
247  	   else
248  	      x.unSetup();
249  	
250  	   ssvec.setSize(0);
251  	   ssvec.forceSetup();
252  	
253  	   solveCount++;
254  	   solveTime->stop();
255  	}
256  	
257  	inline void SLUFactorRational::solveLeft(
258  	   SSVectorRational&      x,
259  	   VectorRational&        y,
260  	   const SVectorRational& rhs1,
261  	   SSVectorRational&      rhs2) //const
262  	{
263  	
264  	   solveTime->start();
265  	
266  	   int   n;
267  	   Rational* svec = ssvec.altValues();
268  	   int*  sidx = ssvec.altIndexMem();
269  	   int   rn   = rhs2.size();
270  	   int*  ridx = rhs2.altIndexMem();
271  	
272  	   x.clear();
273  	   y.clear();
274  	   ssvec.assign(rhs1);
275  	   n = ssvec.size(); // see altValues();
276  	   n = vSolveLeft2(x.altValues(), x.altIndexMem(), svec, sidx, n,
277  	                   y.get_ptr(), rhs2.altValues(), ridx, rn);
278  	
279  	   x.setSize(n);
280  	
281  	   if(n > 0)
282  	      x.forceSetup();
283  	   else
284  	      x.unSetup();
285  	
286  	   rhs2.setSize(0);
287  	   rhs2.forceSetup();
288  	   ssvec.setSize(0);
289  	   ssvec.forceSetup();
290  	
291  	   solveCount++;
292  	   solveTime->stop();
293  	}
294  	
295  	inline void SLUFactorRational::solveLeft(
296  	   SSVectorRational&      x,
297  	   VectorRational&        y,
298  	   VectorRational&        z,
299  	   const SVectorRational& rhs1,
300  	   SSVectorRational&      rhs2,
301  	   SSVectorRational&      rhs3)
302  	{
303  	
304  	   solveTime->start();
305  	
306  	   int   n;
307  	   Rational* svec = ssvec.altValues();
308  	   int*  sidx = ssvec.altIndexMem();
309  	
310  	   x.clear();
311  	   y.clear();
312  	   z.clear();
313  	   ssvec.assign(rhs1);
314  	   n = ssvec.size(); // see altValues();
315  	   n = vSolveLeft3(x.altValues(), x.altIndexMem(), svec, sidx, n,
316  	                   y.get_ptr(), rhs2.altValues(), rhs2.altIndexMem(), rhs2.size(),
317  	                   z.get_ptr(), rhs3.altValues(), rhs3.altIndexMem(), rhs3.size());
318  	
319  	   x.setSize(n);
320  	
321  	   if(n > 0)
322  	      x.forceSetup();
323  	   else
324  	      x.unSetup();
325  	
326  	   ssvec.setSize(0);
327  	   ssvec.forceSetup();
328  	
329  	   solveCount++;
330  	   solveTime->stop();
331  	}
332  	
333  	inline Rational SLUFactorRational::stability() const
334  	{
335  	
336  	   if(status() != OK)
337  	      return 0;
338  	
339  	   if(maxabs < initMaxabs)
340  	      return 1;
341  	
342  	   return initMaxabs / maxabs;
343  	}
344  	
345  	inline std::string SLUFactorRational::statistics() const
346  	{
347  	   std::stringstream s;
348  	   s  << "Factorizations     : " << std::setw(10) << getFactorCount() << std::endl
349  	      << "  Time spent       : " << std::setw(10) << std::fixed << std::setprecision(
350  	         2) << getFactorTime() << std::endl
351  	      << "Solves             : " << std::setw(10) << getSolveCount() << std::endl
352  	      << "  Time spent       : " << std::setw(10) << getSolveTime() << std::endl;
353  	
354  	   return s.str();
355  	}
356  	
357  	inline void SLUFactorRational::changeEta(int idx, SSVectorRational& et)
358  	{
359  	
360  	   int es = et.size(); // see altValues()
361  	   update(idx, et.altValues(), et.altIndexMem(), es);
362  	   et.setSize(0);
363  	   et.forceSetup();
364  	}
365  	
366  	inline SLUFactorRational::Status SLUFactorRational::change(
367  	   int             idx,
368  	   const SVectorRational&  subst,
369  	   const SSVectorRational* e)
370  	{
371  	
372  	   // BH 2005-08-23: The boolean usetup indicates that an "update vector"
373  	   // has been set up. I suppose that SSVectorRational forest is this
374  	   // update vector, which is set up by solveRight4update() and
375  	   // solve2right4update() in order to optimize the basis update.
376  	
377  	   if(usetup)
378  	   {
379  	      if(l.updateType == FOREST_TOMLIN)               // FOREST_TOMLIN updates
380  	      {
381  	         // BH 2005-08-19: The size of a SSVectorRational is the size of the
382  	         // index set, i.e.  the number of nonzeros which is only
383  	         // defined if the SSVectorRational is set up.  Since
384  	         // SSVectorRational::altValues() calls unSetup() the size needs to be
385  	         // stored before the following call.
386  	         int fsize = forest.size(); // see altValues()
387  	         forestUpdate(idx, forest.altValues(), fsize, forest.altIndexMem());
388  	         forest.setSize(0);
389  	         forest.forceSetup();
390  	      }
391  	      else
392  	      {
393  	         assert(l.updateType == ETA);
394  	         changeEta(idx, eta);
395  	      }
396  	   }
397  	   else if(e != 0)                                    // ETA updates
398  	   {
399  	      l.updateType = ETA;
400  	      updateNoClear(idx, e->values(), e->indexMem(), e->size());
401  	      l.updateType = uptype;
402  	   }
403  	   else if(l.updateType == FOREST_TOMLIN)             // FOREST_TOMLIN updates
404  	   {
405  	      assert(0);  // probably this part is never called.
406  	      // forestUpdate() with the last parameter set to NULL should fail.
407  	      forest = subst;
408  	      CLUFactorRational::solveLright(forest.altValues());
409  	      forestUpdate(idx, forest.altValues(), 0, nullptr);
410  	      forest.setSize(0);
411  	      forest.forceSetup();
412  	   }
413  	   else                                               // ETA updates
414  	   {
415  	      assert(l.updateType == ETA);
416  	      vec = subst;
417  	      eta.clear();
418  	      CLUFactorRational::solveRight(eta.altValues(), vec.get_ptr());
419  	      changeEta(idx, eta);
420  	   }
421  	
422  	   usetup = false;
423  	
424  	   SPxOut::debug(this, "DSLUFA01\tupdated\t\tstability = {}\n", stability());
425  	
426  	   return status();
427  	}
428  	
429  	inline void SLUFactorRational::init()
430  	{
431  	   rowMemMult    = 5;          /* factor of minimum Memory * #of nonzeros */
432  	   colMemMult    = 5;          /* factor of minimum Memory * #of nonzeros */
433  	   lMemMult      = 1;          /* factor of minimum Memory * #of nonzeros */
434  	
435  	   l.firstUpdate = 0;
436  	   l.firstUnused = 0;
437  	   thedim        = 0;
438  	
439  	   usetup        = false;
440  	   maxabs        = 1;
441  	   initMaxabs    = 1;
442  	   lastThreshold = minThreshold;
443  	   minStability  = SOPLEX_MINSTABILITYRAT;
444  	   stat          = UNLOADED;
445  	
446  	   vec.clear();
447  	   eta.clear();
448  	   ssvec.clear();
449  	   forest.clear();
450  	
451  	   u.col.size    = 100;
452  	   l.startSize   = 100;
453  	
454  	   l.rval.reDim(0);
455  	
456  	   if(l.ridx)
457  	      spx_free(l.ridx);
458  	
459  	   if(l.rbeg)
460  	      spx_free(l.rbeg);
461  	
462  	   if(l.rorig)
463  	      spx_free(l.rorig);
464  	
465  	   if(l.rperm)
466  	      spx_free(l.rperm);
467  	
468  	   if(u.row.idx)
469  	      spx_free(u.row.idx);
470  	
471  	   if(u.col.idx)
472  	      spx_free(u.col.idx);
473  	
474  	   if(l.idx)
475  	      spx_free(l.idx);
476  	
477  	   if(l.start)
478  	      spx_free(l.start);
479  	
480  	   if(l.row)
481  	      spx_free(l.row);
482  	
483  	   // init() is used in constructor of SLUFactorRational so we have to
484  	   // clean up if anything goes wrong here
485  	   try
486  	   {
487  	      u.row.val.reDim(100);
488  	      spx_alloc(u.row.idx, u.row.val.dim());
489  	      spx_alloc(u.col.idx, u.col.size);
490  	
491  	      l.val.reDim(100);
492  	      spx_alloc(l.idx,   l.val.dim());
493  	      spx_alloc(l.start, l.startSize);
494  	      spx_alloc(l.row,   l.startSize);
495  	   }
496  	   catch(const SPxMemoryException& x)
497  	   {
498  	      freeAll();
499  	      throw x;
500  	   }
501  	}
502  	
503  	inline void SLUFactorRational::clear()
504  	{
505  	   if(stat != UNLOADED)
506  	      init();
507  	}
508  	
509  	/** assignment used to implement operator=() and copy constructor.
510  	 *  If this is initialised, freeAll() has to be called before.
511  	 *  Class objects from SLUFactorRational are not copied here.
512  	 */
513  	inline void SLUFactorRational::assign(const SLUFactorRational& old)
514  	{
515  	   unsigned int thediminc;
516  	
517  	   // slufactor_rational
518  	   uptype        = old.uptype;
519  	   minThreshold  = old.minThreshold;
520  	   minStability  = old.minStability;
521  	   lastThreshold = old.lastThreshold;
522  	
523  	   // clufactor_rational
524  	   stat          = old.stat;
525  	   thedim        = old.thedim;
526  	   nzCnt         = old.nzCnt;
527  	   initMaxabs    = old.initMaxabs;
528  	   maxabs        = old.maxabs;
529  	   rowMemMult    = old.rowMemMult;
530  	   colMemMult    = old.colMemMult;
531  	   lMemMult      = old.lMemMult;
532  	   factorCount   = old.factorCount;
533  	   factorTime    = TimerFactory::createTimer(old.factorTime->type());
534  	   solveTime     = TimerFactory::createTimer(old.solveTime->type());
535  	   timeLimit     = old.timeLimit;
536  	
537  	   spx_alloc(row.perm, thedim);
538  	   spx_alloc(row.orig, thedim);
539  	   spx_alloc(col.perm, thedim);
540  	   spx_alloc(col.orig, thedim);
541  	
542  	   memcpy(row.perm, old.row.perm, (unsigned int)thedim * sizeof(*row.perm));
543  	   memcpy(row.orig, old.row.orig, (unsigned int)thedim * sizeof(*row.orig));
544  	   memcpy(col.perm, old.col.perm, (unsigned int)thedim * sizeof(*col.perm));
545  	   memcpy(col.orig, old.col.orig, (unsigned int)thedim * sizeof(*col.orig));
546  	   diag = old.diag;
547  	
548  	   work = vec.get_ptr();
549  	
550  	   /* setup U
551  	    */
552  	   thediminc = (unsigned int)(thedim + 1);
553  	   u.row.used = old.u.row.used;
554  	   u.row.val = old.u.row.val;
555  	
556  	   spx_alloc(u.row.elem,  thedim);
557  	   spx_alloc(u.row.idx,   u.row.val.dim());
558  	   spx_alloc(u.row.start, thediminc);
559  	   spx_alloc(u.row.len, thediminc);
560  	   spx_alloc(u.row.max, thediminc);
561  	
562  	   memcpy(u.row.elem,  old.u.row.elem, (unsigned int)thedim       * sizeof(*u.row.elem));
563  	   memcpy(u.row.idx,   old.u.row.idx, (unsigned int)u.row.val.dim()   * sizeof(*u.row.idx));
564  	   memcpy(u.row.start, old.u.row.start, thediminc * sizeof(*u.row.start));
565  	   memcpy(u.row.len,   old.u.row.len, thediminc * sizeof(*u.row.len));
566  	   memcpy(u.row.max,   old.u.row.max, thediminc * sizeof(*u.row.max));
567  	
568  	   // need to make row list ok ?
569  	   if(thedim > 0 && stat == OK)
570  	   {
571  	      u.row.list.idx = old.u.row.list.idx; // .idx neu
572  	
573  	      const Dring* oring = &old.u.row.list;
574  	      Dring*       ring  = &u.row.list;
575  	
576  	      while(oring->next != &old.u.row.list)
577  	      {
578  	         ring->next       = &u.row.elem[oring->next->idx];
579  	         ring->next->prev = ring;
580  	         oring            = oring->next;
581  	         ring             = ring->next;
582  	      }
583  	
584  	      ring->next       = &u.row.list;
585  	      ring->next->prev = ring;
586  	   }
587  	
588  	   u.col.size = old.u.col.size;
589  	   u.col.used = old.u.col.used;
590  	
591  	   spx_alloc(u.col.elem,  thedim);
592  	   spx_alloc(u.col.idx,   u.col.size);
593  	   spx_alloc(u.col.start, thediminc);
594  	   spx_alloc(u.col.len, thediminc);
595  	   spx_alloc(u.col.max, thediminc);
596  	   u.col.val = old.u.col.val;
597  	
598  	   memcpy(u.col.elem,  old.u.col.elem, (unsigned int)thedim       * sizeof(*u.col.elem));
599  	   memcpy(u.col.idx,   old.u.col.idx, (unsigned int)u.col.size   * sizeof(*u.col.idx));
600  	   memcpy(u.col.start, old.u.col.start, thediminc * sizeof(*u.col.start));
601  	   memcpy(u.col.len,   old.u.col.len, thediminc * sizeof(*u.col.len));
602  	   memcpy(u.col.max,   old.u.col.max, thediminc * sizeof(*u.col.max));
603  	
604  	   // need to make col list ok
605  	   if(thedim > 0 && stat == OK)
606  	   {
607  	      u.col.list.idx = old.u.col.list.idx; // .idx neu
608  	
609  	      const Dring* oring = &old.u.col.list;
610  	      Dring*       ring  = &u.col.list;
611  	
612  	      while(oring->next != &old.u.col.list)
613  	      {
614  	         ring->next       = &u.col.elem[oring->next->idx];
615  	         ring->next->prev = ring;
616  	         oring            = oring->next;
617  	         ring             = ring->next;
618  	      }
619  	
620  	      ring->next       = &u.col.list;
621  	      ring->next->prev = ring;
622  	   }
623  	
624  	   /* Setup L
625  	    */
626  	   l.startSize   = old.l.startSize;
627  	   l.firstUpdate = old.l.firstUpdate;
628  	   l.firstUnused = old.l.firstUnused;
629  	   l.updateType  = old.l.updateType;
630  	
631  	   l.val = old.l.val;
632  	   spx_alloc(l.idx,   l.val.dim());
633  	   spx_alloc(l.start, l.startSize);
634  	   spx_alloc(l.row,   l.startSize);
635  	
636  	   memcpy(l.idx,   old.l.idx, (unsigned int)l.val.dim() * sizeof(*l.idx));
637  	   memcpy(l.start, old.l.start, (unsigned int)l.startSize * sizeof(*l.start));
638  	   memcpy(l.row,   old.l.row, (unsigned int)l.startSize * sizeof(*l.row));
639  	
640  	   if(old.l.rval.dim() != 0)
641  	   {
642  	      assert(old.l.ridx  != 0);
643  	      assert(old.l.rbeg  != 0);
644  	      assert(old.l.rorig != 0);
645  	      assert(old.l.rperm != 0);
646  	
647  	      int memsize = l.start[l.firstUpdate];
648  	
649  	      l.rval = old.l.rval;
650  	      spx_alloc(l.ridx,  memsize);
651  	      spx_alloc(l.rbeg,  thediminc);
652  	      spx_alloc(l.rorig, thedim);
653  	      spx_alloc(l.rperm, thedim);
654  	
655  	      memcpy(l.ridx,  old.l.ridx, (unsigned int)memsize     * sizeof(*l.ridx));
656  	      memcpy(l.rbeg,  old.l.rbeg, thediminc * sizeof(*l.rbeg));
657  	      memcpy(l.rorig, old.l.rorig, (unsigned int)thedim      * sizeof(*l.rorig));
658  	      memcpy(l.rperm, old.l.rperm, (unsigned int)thedim      * sizeof(*l.rperm));
659  	   }
660  	   else
661  	   {
662  	      assert(old.l.ridx  == 0);
663  	      assert(old.l.rbeg  == 0);
664  	      assert(old.l.rorig == 0);
665  	      assert(old.l.rperm == 0);
666  	
667  	      l.rval.reDim(0);
668  	      l.ridx  = 0;
669  	      l.rbeg  = 0;
670  	      l.rorig = 0;
671  	      l.rperm = 0;
672  	   }
673  	
674  	   assert(row.perm != 0);
675  	   assert(row.orig != 0);
676  	   assert(col.perm != 0);
677  	   assert(col.orig != 0);
678  	
679  	   assert(u.row.elem  != 0);
680  	   assert(u.row.idx   != 0);
681  	   assert(u.row.start != 0);
682  	   assert(u.row.len   != 0);
683  	   assert(u.row.max   != 0);
684  	
685  	   assert(u.col.elem  != 0);
686  	   assert(u.col.idx   != 0);
687  	   assert(u.col.start != 0);
688  	   assert(u.col.len   != 0);
689  	   assert(u.col.max   != 0);
690  	
691  	   assert(l.idx   != 0);
692  	   assert(l.start != 0);
693  	   assert(l.row   != 0);
694  	
695  	}
696  	
697  	inline void SLUFactorRational::freeAll()
698  	{
699  	
700  	   if(row.perm)
701  	      spx_free(row.perm);
702  	
703  	   if(row.orig)
704  	      spx_free(row.orig);
705  	
706  	   if(col.perm)
707  	      spx_free(col.perm);
708  	
709  	   if(col.orig)
710  	      spx_free(col.orig);
711  	
712  	   if(u.row.elem)
713  	      spx_free(u.row.elem);
714  	
715  	   if(u.row.idx)
716  	      spx_free(u.row.idx);
717  	
718  	   if(u.row.start)
719  	      spx_free(u.row.start);
720  	
721  	   if(u.row.len)
722  	      spx_free(u.row.len);
723  	
724  	   if(u.row.max)
725  	      spx_free(u.row.max);
726  	
727  	   if(u.col.elem)
728  	      spx_free(u.col.elem);
729  	
730  	   if(u.col.idx)
731  	      spx_free(u.col.idx);
732  	
733  	   if(u.col.start)
734  	      spx_free(u.col.start);
735  	
736  	   if(u.col.len)
737  	      spx_free(u.col.len);
738  	
739  	   if(u.col.max)
740  	      spx_free(u.col.max);
741  	
742  	   if(l.idx)
743  	      spx_free(l.idx);
744  	
745  	   if(l.start)
746  	      spx_free(l.start);
747  	
748  	   if(l.row)
749  	      spx_free(l.row);
750  	
751  	   if(l.ridx)
752  	      spx_free(l.ridx);
753  	
754  	   if(l.rbeg)
755  	      spx_free(l.rbeg);
756  	
757  	   if(l.rorig)
758  	      spx_free(l.rorig);
759  	
760  	   if(l.rperm)
761  	      spx_free(l.rperm);
762  	
763  	   spx_free(solveTime);
764  	   spx_free(factorTime);
765  	}
766  	
767  	inline SLUFactorRational::~SLUFactorRational()
768  	{
769  	   freeAll();
770  	}
771  	
772  	static Rational betterThreshold(Rational th)
773  	{
774  	   assert(th < 1);
775  	
776  	   if(10 * th < 1)
777  	      th *= 10;
778  	   else if(10 * th < 8)
779  	      th = (th + 1) / 2;
780  	   else if(th < 0.999)
781  	      th = 0.99999;
782  	
783  	   assert(th < 1);
784  	
785  	   return th;
786  	}
787  	
788  	inline SLUFactorRational::Status SLUFactorRational::load(const SVectorRational* matrix[], int dm)
789  	{
790  	   assert(dm     >= 0);
791  	   assert(matrix != 0);
792  	
793  	   Rational lastStability = stability();
794  	
795  	   initDR(u.row.list);
796  	   initDR(u.col.list);
797  	
798  	   usetup        = false;
799  	   l.updateType  = uptype;
800  	   l.firstUpdate = 0;
801  	   l.firstUnused = 0;
802  	
803  	   if(dm != thedim)
804  	   {
805  	      clear();
806  	
807  	      thedim = dm;
808  	      vec.reDim(thedim);
809  	      ssvec.reDim(thedim);
810  	      eta.reDim(thedim);
811  	      forest.reDim(thedim);
812  	      work = vec.get_ptr();
813  	
814  	      spx_realloc(row.perm, thedim);
815  	      spx_realloc(row.orig, thedim);
816  	      spx_realloc(col.perm, thedim);
817  	      spx_realloc(col.orig, thedim);
818  	      diag.reDim(thedim);
819  	
820  	      spx_realloc(u.row.elem,  thedim);
821  	      spx_realloc(u.row.len,   thedim + 1);
822  	      spx_realloc(u.row.max,   thedim + 1);
823  	      spx_realloc(u.row.start, thedim + 1);
824  	
825  	      spx_realloc(u.col.elem,  thedim);
826  	      spx_realloc(u.col.len,   thedim + 1);
827  	      spx_realloc(u.col.max,   thedim + 1);
828  	      spx_realloc(u.col.start, thedim + 1);
829  	
830  	      l.startSize = thedim + SOPLEX_MAXUPDATES;
831  	
832  	      spx_realloc(l.row,   l.startSize);
833  	      spx_realloc(l.start, l.startSize);
834  	   }
835  	   // the last factorization was reasonably stable, so we decrease the Markowitz threshold (stored in lastThreshold) in
836  	   // order favour sparsity
837  	   else if(lastStability > 2.0 * SOPLEX_MINSTABILITYRAT)
838  	   {
839  	      // we reset lastThreshold to its previous value in the sequence minThreshold, betterThreshold(minThreshold),
840  	      // betterThreshold(betterThreshold(minThreshold)), ...
841  	      Rational last   = minThreshold;
842  	      Rational better = betterThreshold(last);
843  	
844  	      while(better < lastThreshold)
845  	      {
846  	         last   = better;
847  	         better = betterThreshold(last);
848  	      }
849  	
850  	      lastThreshold = last;
851  	
852  	      // we reset the minimum stability (which might have been decreased below) to ensure that the increased sparsity
853  	      // does not hurt the stability
854  	      minStability  = 2 * SOPLEX_MINSTABILITYRAT;
855  	   }
856  	
857  	   u.row.list.idx      = thedim;
858  	   u.row.start[thedim] = 0;
859  	   u.row.max[thedim]   = 0;
860  	   u.row.len[thedim]   = 0;
861  	
862  	   u.col.list.idx      = thedim;
863  	   u.col.start[thedim] = 0;
864  	   u.col.max[thedim]   = 0;
865  	   u.col.len[thedim]   = 0;
866  	
867  	   stat = OK;
868  	   factor(matrix, lastThreshold);
869  	
870  	
871  	   SPxOut::debug(this, "DSLUFA02 threshold = {}\tstability = {}\tSOPLEX_MINSTABILITYRAT = {}\n",
872  	                 lastThreshold, stability(), SOPLEX_MINSTABILITYRAT);
873  	   SPX_DEBUG(
874  	      int i;
875  	      FILE* fl = fopen("dump.lp", "w");
876  	      std::cout << "DSLUFA03 Basis:\n";
877  	      int j = 0;
878  	
879  	      for(i = 0; i < dim(); ++i)
880  	      j += matrix[i]->size();
881  	      for(i = 0; i < dim(); ++i)
882  	   {
883  	      for(j = 0; j < matrix[i]->size(); ++j)
884  	            fprintf(fl, "%8d  %8d  %s\n",
885  	                    i + 1, matrix[i]->index(j) + 1, matrix[i]->value(j).str());
886  	      }
887  	   fclose(fl);
888  	   std::cout << "DSLUFA04 LU-Factors:" << std::endl;
889  	             dump();
890  	
891  	             std::cout << "DSLUFA05 threshold = " << lastThreshold
892  	             << "\tstability = " << stability() << std::endl;
893  	   )
894  	
895  	   assert(isConsistent());
896  	   return Status(stat);
897  	}
898  	
899  	
900  	inline bool SLUFactorRational::isConsistent() const
901  	{
902  	#ifdef ENABLE_CONSISTENCY_CHECKS
903  	   return CLUFactorRational::isConsistent();
904  	#else
905  	   return true;
906  	#endif
907  	}
908  	
909  	inline void SLUFactorRational::dump() const
910  	{
911  	   CLUFactorRational::dump();
912  	}
913  	} // namespace soplex
914