Scippy

SoPlex

Sequential object-oriented simPlex

slufactor.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2016 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file slufactor.cpp
17  * @todo SLUfactor seems to be partly an wrapper for CLUFactor (was C).
18  * This should be properly integrated and demangled.
19  * @todo Does is make sense, to call x.clear() when next x.altValues()
20  * is called.
21  */
22 
23 #include <assert.h>
24 #include <sstream>
25 
26 #include "spxdefines.h"
27 #include "slufactor.h"
28 #include "cring.h"
29 #include "spxalloc.h"
30 #include "spxout.h"
31 #include "exceptions.h"
32 
33 #ifdef SOPLEX_DEBUG
34 #include <stdio.h>
35 #endif
36 
37 namespace soplex
38 {
39 #define MINSTABILITY REAL(4e-2)
40 
41 void SLUFactor::solveRight(Vector& x, const Vector& b) //const
42 {
43 
44  solveTime->start();
45 
46  vec = b;
48 
49  solveCount++;
50  solveTime->stop();
51 }
52 
53 void SLUFactor::solveRight(SSVector& x, const SVector& b) //const
54 {
55 
56  solveTime->start();
57 
58  vec.assign(b);
59  x.clear();
61 
62  solveCount++;
63  solveTime->stop();
64 }
65 
67 {
68 
69  solveTime->start();
70 
71  int m;
72  int n;
73  int f;
74 
75  x.clear();
76  ssvec = b;
77  n = ssvec.size();
78  if (l.updateType == ETA)
79  {
81  ssvec.altValues(), ssvec.altIndexMem(), n, 0, 0, 0);
82  x.setSize(m);
83  //x.forceSetup();
84  x.unSetup();
86  }
87  else
88  {
89  forest.clear();
93  forest.setSize(f);
95  x.setSize(m);
96  x.forceSetup();
97  }
98  usetup = true;
99  ssvec.setSize(0);
100  ssvec.forceSetup();
101 
102  solveCount++;
103  solveTime->stop();
104 }
105 
107  SSVector& x,
108  Vector& y,
109  const SVector& b,
110  SSVector& rhs)
111 {
112 
113  solveTime->start();
114 
115  int m;
116  int n;
117  int f;
118  int* sidx = ssvec.altIndexMem();
119  ssvec.setSize(0);
120  ssvec.forceSetup();
121  int rsize = rhs.size();
122  int* ridx = rhs.altIndexMem();
123 
124  x.clear();
125  y.clear();
126  usetup = true;
127  ssvec = b;
128 
129  if (l.updateType == ETA)
130  {
131  n = ssvec.size();
133  ssvec.get_ptr(), sidx, n, y.get_ptr(),
134  rhs.getEpsilon(), rhs.altValues(), ridx, rsize, 0, 0, 0);
135  x.setSize(m);
136  // x.forceSetup();
137  x.unSetup();
139  }
140  else
141  {
142  forest.clear();
143  n = ssvec.size();
145  ssvec.get_ptr(), sidx, n, y.get_ptr(),
146  rhs.getEpsilon(), rhs.altValues(), ridx, rsize,
148  x.setSize(m);
149  x.forceSetup();
150  forest.setSize(f);
151  forest.forceSetup();
152  }
153  rhs.forceSetup();
154  ssvec.setSize(0);
155  ssvec.forceSetup();
156 
157  solveCount += 2;
158  solveTime->stop();
159 }
160 
162  SSVector& x,
163  SSVector& y,
164  const SVector& b,
165  SSVector& rhs)
166 {
167 
168  solveTime->start();
169 
170  int n;
171  int f;
172  int* sidx = ssvec.altIndexMem();
173  ssvec.setSize(0);
174  ssvec.forceSetup();
175  int rsize = rhs.size();
176  int* ridx = rhs.altIndexMem();
177 
178  x.clear();
179  y.clear();
180  usetup = true;
181  ssvec = b;
182 
183  if (l.updateType == ETA)
184  {
185  n = ssvec.size();
187  ssvec.get_ptr(), sidx, n,
188  y.getEpsilon(), y.altValues(), y.altIndexMem(),
189  rhs.altValues(), ridx, rsize,
190  0, 0, 0);
191  x.setSize(n);
192  // x.forceSetup();
193  x.unSetup();
194  y.setSize(rsize);
195  y.unSetup();
197  }
198  else
199  {
200  forest.clear();
201  n = ssvec.size();
203  ssvec.get_ptr(), sidx, n,
204  y.getEpsilon(), y.altValues(), y.altIndexMem(),
205  rhs.altValues(), ridx, rsize,
207  x.setSize(n);
208  x.forceSetup();
209  y.setSize(rsize);
210  y.forceSetup();
211  forest.setSize(f);
212  forest.forceSetup();
213  }
214  rhs.forceSetup();
215  ssvec.setSize(0);
216  ssvec.forceSetup();
217 
218  solveCount += 2;
219  solveTime->stop();
220 }
221 
222 
224  SSVector& x,
225  Vector& y,
226  Vector& y2,
227  const SVector& b,
228  SSVector& rhs,
229  SSVector& rhs2)
230 {
231 
232  solveTime->start();
233 
234  int m;
235  int n;
236  int f;
237  int* sidx = ssvec.altIndexMem();
238  ssvec.setSize(0);
239  ssvec.forceSetup();
240  int rsize = rhs.size();
241  int* ridx = rhs.altIndexMem();
242  int rsize2 = rhs2.size();
243  int* ridx2 = rhs2.altIndexMem();
244 
245  x.clear();
246  y.clear();
247  y2.clear();
248  usetup = true;
249  ssvec = b;
250 
251  if (l.updateType == ETA)
252  {
253  n = ssvec.size();
255  x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n,
256  y.get_ptr(), rhs.getEpsilon(), rhs.altValues(), ridx, rsize,
257  y2.get_ptr(), rhs2.getEpsilon(),rhs2.altValues(), ridx2, rsize2,
258  0, 0, 0);
259  x.setSize(m);
260  // x.forceSetup();
261  x.unSetup();
263  }
264  else
265  {
266  forest.clear();
267  n = ssvec.size();
269  x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n,
270  y.get_ptr(), rhs.getEpsilon(), rhs.altValues(), ridx, rsize,
271  y2.get_ptr(), rhs2.getEpsilon(),rhs2.altValues(), ridx2, rsize2,
273  x.setSize(m);
274  x.forceSetup();
275  forest.setSize(f);
276  forest.forceSetup();
277  }
278  rhs.forceSetup();
279  rhs2.forceSetup();
280  ssvec.setSize(0);
281  ssvec.forceSetup();
282 
283  solveCount += 3;
284  solveTime->stop();
285 }
286 
288  SSVector& x,
289  SSVector& y,
290  SSVector& y2,
291  const SVector& b,
292  SSVector& rhs,
293  SSVector& rhs2)
294 {
295  solveTime->start();
296 
297  int n;
298  int f;
299  int* sidx = ssvec.altIndexMem();
300  ssvec.setSize(0);
301  ssvec.forceSetup();
302  int rsize = rhs.size();
303  int* ridx = rhs.altIndexMem();
304  int rsize2 = rhs2.size();
305  int* ridx2 = rhs2.altIndexMem();
306 
307  x.clear();
308  y.clear();
309  y2.clear();
310  usetup = true;
311  ssvec = b;
312 
313  if (l.updateType == ETA)
314  {
315  n = ssvec.size();
317  ssvec.get_ptr(), sidx, n,
318  y.getEpsilon(), y.altValues(), y.altIndexMem(),
319  rhs.altValues(), ridx, rsize,
320  y2.getEpsilon(), y2.altValues(), y2.altIndexMem(),
321  rhs2.altValues(), ridx2, rsize2,
322  0, 0, 0);
323  x.setSize(n);
324  // x.forceSetup();
325  x.unSetup();
326  y.setSize(rsize);
327  y.unSetup();
328  y2.setSize(rsize2);
329  y2.unSetup();
331  }
332  else
333  {
334  forest.clear();
335  n = ssvec.size();
337  ssvec.get_ptr(), sidx, n,
338  y.getEpsilon(), y.altValues(), y.altIndexMem(),
339  rhs.altValues(), ridx, rsize,
340  y2.getEpsilon(), y2.altValues(), y2.altIndexMem(),
341  rhs2.altValues(), ridx2, rsize2,
343  x.setSize(n);
344  x.forceSetup();
345  y.setSize(rsize);
346  y.forceSetup();
347  y2.setSize(rsize2);
348  y2.forceSetup();
349 
350  forest.setSize(f);
351  forest.forceSetup();
352  }
353  rhs.forceSetup();
354  rhs2.forceSetup();
355  ssvec.setSize(0);
356  ssvec.forceSetup();
357 
358  solveCount += 3;
359  solveTime->stop();
360 }
361 
362 
363 void SLUFactor::solveLeft(Vector& x, const Vector& b) //const
364 {
365 
366  solveTime->start();
367 
368  vec = b;
369  ///@todo Why is x.clear() here used and not with solveRight() ?
370  x.clear();
372 
373  solveCount++;
374  solveTime->stop();
375 }
376 
377 void SLUFactor::solveLeft(SSVector& x, const SVector& b) //const
378 {
379 
380  solveTime->start();
381 
382  // copy to SSVec is done to avoid having to deal with the Nonzero datatype
383  // TODO change SVec to standard sparse format
384  ssvec.assign(b);
385 
386  x.clear();
387  int sz = ssvec.size(); // see .altValues()
388  int n = vSolveLeft(x.getEpsilon(), x.altValues(), x.altIndexMem(),
389  ssvec.altValues(), ssvec.altIndexMem(), sz);
390 
391  if (n > 0)
392  {
393  x.setSize(n);
394  x.forceSetup();
395  }
396  else
397  x.unSetup();
398 
399  ssvec.setSize(0);
400  ssvec.forceSetup();
401 
402  solveCount++;
403  solveTime->stop();
404 }
405 
407  SSVector& x,
408  Vector& y,
409  const SVector& rhs1,
410  SSVector& rhs2) //const
411 {
412 
413  solveTime->start();
414 
415  int n;
416  Real* svec = ssvec.altValues();
417  int* sidx = ssvec.altIndexMem();
418  int rn = rhs2.size();
419  int* ridx = rhs2.altIndexMem();
420 
421  x.clear();
422  y.clear();
423  ssvec.assign(rhs1);
424  n = ssvec.size(); // see altValues();
425  n = vSolveLeft2(x.getEpsilon(), x.altValues(), x.altIndexMem(), svec, sidx, n,
426  y.get_ptr(), rhs2.altValues(), ridx, rn);
427 
428  // this will unsetup x
429  x.setSize(n);
430 
431  if (n > 0)
432  x.forceSetup();
433 
434  ssvec.setSize(0);
435  ssvec.forceSetup();
436 
437  solveCount += 2;
438  solveTime->stop();
439 }
440 
442  SSVector& x,
443  SSVector& y,
444  const SVector& rhs1,
445  SSVector& rhs2) //const
446 {
447 
448  solveTime->start();
449 
450  int n1, n2;
451  Real* svec = ssvec.altValues();
452  int* sidx = ssvec.altIndexMem();
453 
454  x.clear();
455  y.clear();
456  ssvec.assign(rhs1);
457  n1 = ssvec.size(); // see altValues();
458  n2 = rhs2.size();
459  if( n2 < 10 )
460  {
462  x.altValues(), x.altIndexMem(),
463  svec, sidx, n1,
464  y.altValues(), y.altIndexMem(),
465  rhs2.altValues(), rhs2.altIndexMem(), n2);
466  y.setSize(n2);
467  if( n2 > 0 )
468  y.forceSetup();
469  }
470  else
471  {
472  n1 = vSolveLeft2(x.getEpsilon(), x.altValues(), x.altIndexMem(), svec, sidx, n1,
473  y.altValues(), rhs2.altValues(), rhs2.altIndexMem(), n2);
474 // y.setup();
475  }
476  x.setSize(n1);
477 
478  if (n1 > 0)
479  x.forceSetup();
480 // if (n2 > 0)
481 // y.forceSetup();
482 
483  ssvec.setSize(0);
484  ssvec.forceSetup();
485 
486  solveCount += 2;
487  solveTime->stop();
488 }
489 
490 
492  SSVector& x,
493  Vector& y,
494  Vector& z,
495  const SVector& rhs1,
496  SSVector& rhs2,
497  SSVector& rhs3)
498 {
499 
500  solveTime->start();
501 
502  int n, n2, n3;
503  Real* svec = ssvec.altValues();
504  int* sidx = ssvec.altIndexMem();
505 
506  x.clear();
507  y.clear();
508  z.clear();
509  ssvec.assign(rhs1);
510  n = ssvec.size(); // see altValues();
511  n2 = rhs2.size();
512  n3 = rhs3.size();
513 
514  n = vSolveLeft3(x.getEpsilon(), x.altValues(), x.altIndexMem(), svec, sidx, n,
515  y.get_ptr(), rhs2.altValues(), rhs2.altIndexMem(), n2,
516  z.get_ptr(), rhs3.altValues(), rhs3.altIndexMem(), n3);
517 
518  x.setSize(n);
519 
520  if (n > 0)
521  x.forceSetup();
522 
523  ssvec.setSize(0);
524  ssvec.forceSetup();
525 
526  solveCount += 3;
527  solveTime->stop();
528 }
529 
531  SSVector& x,
532  SSVector& y,
533  SSVector& z,
534  const SVector& rhs1,
535  SSVector& rhs2,
536  SSVector& rhs3)
537 {
538 
539  solveTime->start();
540 
541  int n1, n2, n3;
542  Real* svec = ssvec.altValues();
543  int* sidx = ssvec.altIndexMem();
544 
545  x.clear();
546  y.clear();
547  z.clear();
548  ssvec.assign(rhs1);
549  n1 = ssvec.size(); // see altValues();
550  n2 = rhs2.size();
551  n3 = rhs3.size();
553  x.altValues(), x.altIndexMem(),
554  svec, sidx, n1,
555  y.altValues(), y.altIndexMem(),
556  rhs2.altValues(), rhs2.altIndexMem(), n2,
557  z.altValues(), z.altIndexMem(),
558  rhs3.altValues(), rhs3.altIndexMem(), n3);
559  x.setSize(n1);
560  y.setSize(n2);
561  z.setSize(n3);
562 
563  if (n1 > 0)
564  x.forceSetup();
565  if (n2 > 0)
566  y.forceSetup();
567  if (n3 > 0)
568  z.forceSetup();
569 
570  ssvec.setSize(0);
571  ssvec.forceSetup();
572 
573  solveCount += 3;
574  solveTime->stop();
575 }
576 
577 
579 {
580 
581  if (status() != OK)
582  return 0;
583 
584  if (maxabs < initMaxabs)
585  return 1;
586 
587  return initMaxabs / maxabs;
588 }
589 
590 std::string SLUFactor::statistics() const
591 {
592  std::stringstream s;
593  s << "Factorizations : " << std::setw(10) << getFactorCount() << std::endl
594  << " Time spent : " << std::setw(10) << std::fixed << std::setprecision(2) << getFactorTime() << std::endl
595  << "Solves : " << std::setw(10) << getSolveCount() << std::endl
596  << " Time spent : " << std::setw(10) << getSolveTime() << std::endl;
597 
598  return s.str();
599 }
600 
601 void SLUFactor::changeEta(int idx, SSVector& et)
602 {
603 
604  int es = et.size(); // see altValues()
605  update(idx, et.altValues(), et.altIndexMem(), es);
606  et.setSize(0);
607  et.forceSetup();
608 }
609 
611  int idx,
612  const SVector& subst,
613  const SSVector* e)
614 {
615 
616  // BH 2005-08-23: The boolean usetup indicates that an "update vector"
617  // has been set up. I suppose that SSVector forest is this
618  // update vector, which is set up by solveRight4update() and
619  // solve2right4update() in order to optimize the basis update.
620 
621  if (usetup)
622  {
623  if (l.updateType == FOREST_TOMLIN) // FOREST_TOMLIN updates
624  {
625  // BH 2005-08-19: The size of a SSVector is the size of the
626  // index set, i.e. the number of nonzeros which is only
627  // defined if the SSVector is set up. Since
628  // SSVector::altValues() calls unSetup() the size needs to be
629  // stored before the following call.
630  int fsize = forest.size(); // see altValues()
631  forestUpdate(idx, forest.altValues(), fsize, forest.altIndexMem());
632  forest.setSize(0);
633  forest.forceSetup();
634  }
635  else
636  {
637  assert(l.updateType == ETA);
638  changeEta(idx, eta);
639  }
640  }
641  else if (e != 0) // ETA updates
642  {
643  l.updateType = ETA;
644  updateNoClear(idx, e->values(), e->indexMem(), e->size());
645  l.updateType = uptype;
646  }
647  else if (l.updateType == FOREST_TOMLIN) // FOREST_TOMLIN updates
648  {
649  assert(0); // probably this part is never called.
650  // forestUpdate() with the last parameter set to NULL should fail.
651  forest = subst;
653  forestUpdate(idx, forest.altValues(), 0, 0);
654  forest.setSize(0);
655  forest.forceSetup();
656  }
657  else // ETA updates
658  {
659  assert(l.updateType == ETA);
660  vec = subst;
661  eta.clear();
663  changeEta(idx, eta);
664  }
665  usetup = false;
666 
667  MSG_DEBUG( std::cout << "DSLUFA01\tupdated\t\tstability = " << stability()
668  << std::endl; )
669 
670  return status();
671 }
672 
674 {
675 
676  rowMemMult = 5; /* factor of minimum Memory * #of nonzeros */
677  colMemMult = 5; /* factor of minimum Memory * #of nonzeros */
678  lMemMult = 1; /* factor of minimum Memory * #of nonzeros */
679 
680  l.firstUpdate = 0;
681  l.firstUnused = 0;
682  thedim = 0;
683 
685  usetup = false;
686  maxabs = 1;
687  initMaxabs = 1;
690  stat = UNLOADED;
691 
692  vec.clear();
693  eta.clear();
694  ssvec.clear();
695  forest.clear();
696 
697  u.row.size = 100;
698  u.col.size = 100;
699  l.size = 100;
700  l.startSize = 100;
701 
702  if (l.rval)
703  spx_free(l.rval);
704  if(l.ridx)
705  spx_free(l.ridx);
706  if(l.rbeg)
707  spx_free(l.rbeg);
708  if(l.rorig)
709  spx_free(l.rorig);
710  if(l.rperm)
711  spx_free(l.rperm);
712 
713  if(u.row.val)
714  spx_free(u.row.val);
715  if(u.row.idx)
716  spx_free(u.row.idx);
717  if(u.col.idx)
718  spx_free(u.col.idx);
719  if(l.val)
720  spx_free(l.val);
721  if(l.idx)
722  spx_free(l.idx);
723  if(l.start)
724  spx_free(l.start);
725  if(l.row)
726  spx_free(l.row);
727 
728  // G clear() is used in constructor of SLUFactor so we have to
729  // G clean up if anything goes wrong here
730  try
731  {
732  spx_alloc(u.row.val, u.row.size);
733  spx_alloc(u.row.idx, u.row.size);
734  spx_alloc(u.col.idx, u.col.size);
735 
736  spx_alloc(l.val, l.size);
737  spx_alloc(l.idx, l.size);
740  }
741  catch(const SPxMemoryException& x)
742  {
743  freeAll();
744  throw x;
745  }
746 }
747 
748 /** assignment used to implement operator=() and copy constructor.
749  * If this is initialised, freeAll() has to be called before.
750  * Class objects from SLUFactor are not copied here.
751  */
752 void SLUFactor::assign(const SLUFactor& old)
753 {
754  spxout = old.spxout;
755 
756  // slufactor
757  uptype = old.uptype;
760  epsilon = old.epsilon;
762 
763  // clufactor
764  stat = old.stat;
765  thedim = old.thedim;
766  nzCnt = old.nzCnt;
767  initMaxabs = old.initMaxabs;
768  maxabs = old.maxabs;
769  rowMemMult = old.rowMemMult;
770  colMemMult = old.colMemMult;
771  lMemMult = old.lMemMult;
772 
778 
779  memcpy(row.perm, old.row.perm, (unsigned int)thedim * sizeof(*row.perm));
780  memcpy(row.orig, old.row.orig, (unsigned int)thedim * sizeof(*row.orig));
781  memcpy(col.perm, old.col.perm, (unsigned int)thedim * sizeof(*col.perm));
782  memcpy(col.orig, old.col.orig, (unsigned int)thedim * sizeof(*col.orig));
783  memcpy(diag, old.diag, (unsigned int)thedim * sizeof(*diag));
784 
785  work = vec.get_ptr();
786 
787  /* setup U
788  */
789  u.row.size = old.u.row.size;
790  u.row.used = old.u.row.used;
791 
793  spx_alloc(u.row.val, u.row.size);
794  spx_alloc(u.row.idx, u.row.size);
795  spx_alloc(u.row.start, thedim + 1);
796  spx_alloc(u.row.len, thedim + 1);
797  spx_alloc(u.row.max, thedim + 1);
798 
799  memcpy(u.row.elem, old.u.row.elem, (unsigned int)thedim * sizeof(*u.row.elem));
800  memcpy(u.row.val, old.u.row.val, (unsigned int)u.row.size * sizeof(*u.row.val));
801  memcpy(u.row.idx, old.u.row.idx, (unsigned int)u.row.size * sizeof(*u.row.idx));
802  memcpy(u.row.start, old.u.row.start, (unsigned int)(thedim + 1) * sizeof(*u.row.start));
803  memcpy(u.row.len, old.u.row.len, (unsigned int)(thedim + 1) * sizeof(*u.row.len));
804  memcpy(u.row.max, old.u.row.max, (unsigned int)(thedim + 1) * sizeof(*u.row.max));
805 
806  // need to make row list ok ?
807  if (thedim > 0 && stat == OK)
808  {
809  u.row.list.idx = old.u.row.list.idx; // .idx neu
810 
811  const Dring* oring = &old.u.row.list;
812  Dring* ring = &u.row.list;
813 
814  while(oring->next != &old.u.row.list)
815  {
816  ring->next = &u.row.elem[oring->next->idx];
817  ring->next->prev = ring;
818  oring = oring->next;
819  ring = ring->next;
820  }
821  ring->next = &u.row.list;
822  ring->next->prev = ring;
823  }
824 
825  u.col.size = old.u.col.size;
826  u.col.used = old.u.col.used;
827 
829  spx_alloc(u.col.idx, u.col.size);
830  spx_alloc(u.col.start, thedim + 1);
831  spx_alloc(u.col.len, thedim + 1);
832  spx_alloc(u.col.max, thedim + 1);
833 
834  if (old.u.col.val != 0)
835  {
836  spx_alloc(u.col.val, u.col.size);
837  memcpy(u.col.val, old.u.col.val, (unsigned int)u.col.size * sizeof(*u.col.val));
838  }
839  else
840  u.col.val = 0;
841 
842  memcpy(u.col.elem, old.u.col.elem, (unsigned int)thedim * sizeof(*u.col.elem));
843  memcpy(u.col.idx, old.u.col.idx, (unsigned int)u.col.size * sizeof(*u.col.idx));
844  memcpy(u.col.start, old.u.col.start, (unsigned int)(thedim + 1) * sizeof(*u.col.start));
845  memcpy(u.col.len, old.u.col.len, (unsigned int)(thedim + 1) * sizeof(*u.col.len));
846  memcpy(u.col.max, old.u.col.max, (unsigned int)(thedim + 1) * sizeof(*u.col.max));
847 
848  // need to make col list ok
849  if (thedim > 0 && stat == OK)
850  {
851  u.col.list.idx = old.u.col.list.idx; // .idx neu
852 
853  const Dring* oring = &old.u.col.list;
854  Dring* ring = &u.col.list;
855 
856  while(oring->next != &old.u.col.list)
857  {
858  ring->next = &u.col.elem[oring->next->idx];
859  ring->next->prev = ring;
860  oring = oring->next;
861  ring = ring->next;
862  }
863  ring->next = &u.col.list;
864  ring->next->prev = ring;
865  }
866 
867  /* Setup L
868  */
869  l.size = old.l.size;
870  l.startSize = old.l.startSize;
871  l.firstUpdate = old.l.firstUpdate;
872  l.firstUnused = old.l.firstUnused;
873  l.updateType = old.l.updateType;
874 
875  spx_alloc(l.val, l.size);
876  spx_alloc(l.idx, l.size);
879 
880  memcpy(l.val, old.l.val, (unsigned int)l.size * sizeof(*l.val));
881  memcpy(l.idx, old.l.idx, (unsigned int)l.size * sizeof(*l.idx));
882  memcpy(l.start, old.l.start, (unsigned int)l.startSize * sizeof(*l.start));
883  memcpy(l.row, old.l.row, (unsigned int)l.startSize * sizeof(*l.row));
884 
885  if (old.l.rval != 0)
886  {
887  assert(old.l.ridx != 0);
888  assert(old.l.rbeg != 0);
889  assert(old.l.rorig != 0);
890  assert(old.l.rperm != 0);
891 
892  int memsize = l.start[l.firstUpdate];
893 
894  spx_alloc(l.rval, memsize);
895  spx_alloc(l.ridx, memsize);
896  spx_alloc(l.rbeg, thedim + 1);
899 
900  memcpy(l.rval, old.l.rval, (unsigned int)memsize * sizeof(*l.rval));
901  memcpy(l.ridx, old.l.ridx, (unsigned int)memsize * sizeof(*l.ridx));
902  memcpy(l.rbeg, old.l.rbeg, (unsigned int)(thedim + 1) * sizeof(*l.rbeg));
903  memcpy(l.rorig, old.l.rorig, (unsigned int)thedim * sizeof(*l.rorig));
904  memcpy(l.rperm, old.l.rperm, (unsigned int)thedim * sizeof(*l.rperm));
905  }
906  else
907  {
908  assert(old.l.ridx == 0);
909  assert(old.l.rbeg == 0);
910  assert(old.l.rorig == 0);
911  assert(old.l.rperm == 0);
912 
913  l.rval = 0;
914  l.ridx = 0;
915  l.rbeg = 0;
916  l.rorig = 0;
917  l.rperm = 0;
918  }
919 
920  assert(row.perm != 0);
921  assert(row.orig != 0);
922  assert(col.perm != 0);
923  assert(col.orig != 0);
924  assert(diag != 0);
925 
926  assert(u.row.elem != 0);
927  assert(u.row.val != 0);
928  assert(u.row.idx != 0);
929  assert(u.row.start != 0);
930  assert(u.row.len != 0);
931  assert(u.row.max != 0);
932 
933  assert(u.col.elem != 0);
934  assert(u.col.idx != 0);
935  assert(u.col.start != 0);
936  assert(u.col.len != 0);
937  assert(u.col.max != 0);
938 
939  assert(l.val != 0);
940  assert(l.idx != 0);
941  assert(l.start != 0);
942  assert(l.row != 0);
943 
944 }
945 
947 {
948 
949  if (this != &old)
950  {
951  // we don't need to copy them, because they are temporary vectors
952  vec.clear();
953  ssvec.clear();
954 
955  eta = old.eta;
956  forest = old.forest;
957 
958  timerType = old.timerType;
959 
960  freeAll();
961  try
962  {
963  assign(old);
964  }
965  catch(const SPxMemoryException& x)
966  {
967  freeAll();
968  throw x;
969  }
970  assert(isConsistent());
971  }
972  return *this;
973 }
974 
976  : CLUFactor()
977  , vec (1)
978  , ssvec (1)
979  , usetup (false)
981  , eta (1)
982  , forest (1)
983  , minThreshold (0.01)
984  , timerType(Timer::USER_TIME)
985 {
986  row.perm = 0;
987  row.orig = 0;
988  col.perm = 0;
989  col.orig = 0;
990  diag = 0;
991  u.row.elem = 0;
992  u.row.val = 0;
993  u.row.idx = 0;
994  u.row.start = 0;
995  u.row.len = 0;
996  u.row.max = 0;
997  u.col.elem = 0;
998  u.col.idx = 0;
999  u.col.start = 0;
1000  u.col.len = 0;
1001  u.col.max = 0;
1002  u.col.val = 0;
1003  l.val = 0;
1004  l.idx = 0;
1005  l.start = 0;
1006  l.row = 0;
1007  l.rval = 0;
1008  l.ridx = 0;
1009  l.rbeg = 0;
1010  l.rorig = 0;
1011  l.rperm = 0;
1012 
1013  nzCnt = 0;
1014  thedim = 0;
1015  try
1016  {
1023  spx_alloc(diag, thedim);
1024 
1025  work = vec.get_ptr();
1026 
1027  u.row.size = 1;
1028  u.row.used = 0;
1030  spx_alloc(u.row.val, u.row.size);
1031  spx_alloc(u.row.idx, u.row.size);
1032  spx_alloc(u.row.start, thedim + 1);
1033  spx_alloc(u.row.len, thedim + 1);
1034  spx_alloc(u.row.max, thedim + 1);
1035 
1036  u.row.list.idx = thedim;
1037  u.row.start[thedim] = 0;
1038  u.row.max [thedim] = 0;
1039  u.row.len [thedim] = 0;
1040 
1041  u.col.size = 1;
1042  u.col.used = 0;
1044  spx_alloc(u.col.idx, u.col.size);
1045  spx_alloc(u.col.start, thedim + 1);
1046  spx_alloc(u.col.len, thedim + 1);
1047  spx_alloc(u.col.max, thedim + 1);
1048  u.col.val = 0;
1049 
1050  u.col.list.idx = thedim;
1051  u.col.start[thedim] = 0;
1052  u.col.max[thedim] = 0;
1053  u.col.len[thedim] = 0;
1054 
1055  l.size = 1;
1056 
1057  spx_alloc(l.val, l.size);
1058  spx_alloc(l.idx, l.size);
1059 
1060  l.startSize = 1;
1061  l.firstUpdate = 0;
1062  l.firstUnused = 0;
1063 
1066  }
1067  catch(const SPxMemoryException& x)
1068  {
1069  freeAll();
1070  throw x;
1071  }
1072 
1073  l.rval = 0;
1074  l.ridx = 0;
1075  l.rbeg = 0;
1076  l.rorig = 0;
1077  l.rperm = 0;
1078 
1079  SLUFactor::clear(); // clear() is virtual
1080 
1081  factorCount = 0;
1082  solveCount = 0;
1083 
1084  assert(row.perm != 0);
1085  assert(row.orig != 0);
1086  assert(col.perm != 0);
1087  assert(col.orig != 0);
1088  assert(diag != 0);
1089 
1090  assert(u.row.elem != 0);
1091  assert(u.row.val != 0);
1092  assert(u.row.idx != 0);
1093  assert(u.row.start != 0);
1094  assert(u.row.len != 0);
1095  assert(u.row.max != 0);
1096 
1097  assert(u.col.elem != 0);
1098  assert(u.col.idx != 0);
1099  assert(u.col.start != 0);
1100  assert(u.col.len != 0);
1101  assert(u.col.max != 0);
1102 
1103  assert(l.val != 0);
1104  assert(l.idx != 0);
1105  assert(l.start != 0);
1106  assert(l.row != 0);
1107 
1108  assert(SLUFactor::isConsistent());
1109 }
1110 
1112  : SLinSolver( old )
1113  , CLUFactor()
1114  , vec(1) // we don't need to copy it, because they are temporary vectors
1115  , ssvec(1) // we don't need to copy it, because they are temporary vectors
1116  , usetup(old.usetup)
1117  , eta (old.eta)
1118  , forest(old.forest)
1119  , timerType(old.timerType)
1120 {
1121  row.perm = 0;
1122  row.orig = 0;
1123  col.perm = 0;
1124  col.orig = 0;
1125  diag = 0;
1126  u.row.elem = 0;
1127  u.row.val = 0;
1128  u.row.idx = 0;
1129  u.row.start = 0;
1130  u.row.len = 0;
1131  u.row.max = 0;
1132  u.col.elem = 0;
1133  u.col.idx = 0;
1134  u.col.start = 0;
1135  u.col.len = 0;
1136  u.col.max = 0;
1137  u.col.val = 0;
1138  l.val = 0;
1139  l.idx = 0;
1140  l.start = 0;
1141  l.row = 0;
1142  l.rval = 0;
1143  l.ridx = 0;
1144  l.rbeg = 0;
1145  l.rorig = 0;
1146  l.rperm = 0;
1147 
1148  solveCount = 0;
1151 
1152  try
1153  {
1154  assign(old);
1155  }
1156  catch(const SPxMemoryException& x)
1157  {
1158  freeAll();
1159  throw x;
1160  }
1161  assert(SLUFactor::isConsistent());
1162 }
1163 
1165 {
1166 
1167  if(row.perm) spx_free(row.perm);
1168  if(row.orig) spx_free(row.orig);
1169  if(col.perm) spx_free(col.perm);
1170  if(col.orig) spx_free(col.orig);
1171  if(u.row.elem) spx_free(u.row.elem);
1172  if(u.row.val) spx_free(u.row.val);
1173  if(u.row.idx) spx_free(u.row.idx);
1174  if(u.row.start) spx_free(u.row.start);
1175  if(u.row.len) spx_free(u.row.len);
1176  if(u.row.max) spx_free(u.row.max);
1177  if(u.col.elem) spx_free(u.col.elem);
1178  if(u.col.idx) spx_free(u.col.idx);
1179  if(u.col.start) spx_free(u.col.start);
1180  if(u.col.len) spx_free(u.col.len);
1181  if(u.col.max) spx_free(u.col.max);
1182  if(l.val) spx_free(l.val);
1183  if(l.idx) spx_free(l.idx);
1184  if(l.start) spx_free(l.start);
1185  if(l.row) spx_free(l.row);
1186 
1187  if(diag) spx_free(diag);
1188 
1189  if (u.col.val) spx_free(u.col.val);
1190 
1191  if (l.rval) spx_free(l.rval);
1192  if(l.ridx) spx_free(l.ridx);
1193  if(l.rbeg) spx_free(l.rbeg);
1194  if(l.rorig) spx_free(l.rorig);
1195  if(l.rperm) spx_free(l.rperm);
1196 }
1197 
1199 {
1200  freeAll();
1201 
1202  solveTime->~Timer();
1203  factorTime->~Timer();
1206 }
1207 
1209 {
1210  assert(th < REAL(1.0));
1211 
1212  if (LT(th, REAL(0.1)))
1213  th *= REAL(10.0);
1214  else if (LT(th, REAL(0.9)))
1215  th = (th + REAL(1.0)) / REAL(2.0);
1216  else if (LT(th, REAL(0.999)))
1217  th = REAL(0.99999);
1218 
1219  assert(th < REAL(1.0));
1220 
1221  return th;
1222 }
1223 
1225 {
1226  assert(dm >= 0);
1227  assert(matrix != 0);
1228 
1229  Real lastStability = stability();
1230 
1231  initDR(u.row.list);
1232  initDR(u.col.list);
1233 
1234  usetup = false;
1235  l.updateType = uptype;
1236  l.firstUpdate = 0;
1237  l.firstUnused = 0;
1238 
1239  if (dm != thedim)
1240  {
1241  clear();
1242 
1243  thedim = dm;
1244  vec.reDim(thedim);
1245  ssvec.reDim(thedim);
1246  eta.reDim(thedim);
1247  forest.reDim(thedim);
1248  work = vec.get_ptr();
1249 
1255 
1257  spx_realloc(u.row.len, thedim + 1);
1258  spx_realloc(u.row.max, thedim + 1);
1259  spx_realloc(u.row.start, thedim + 1);
1260 
1262  spx_realloc(u.col.len, thedim + 1);
1263  spx_realloc(u.col.max, thedim + 1);
1264  spx_realloc(u.col.start, thedim + 1);
1265 
1267 
1270  }
1271  // the last factorization was reasonably stable, so we decrease the Markowitz threshold (stored in lastThreshold) in
1272  // order favour sparsity
1273  else if (lastStability > 2.0 * minStability)
1274  {
1275  // we reset lastThreshold to its previous value in the sequence minThreshold, betterThreshold(minThreshold),
1276  // betterThreshold(betterThreshold(minThreshold)), ...
1277  Real last = minThreshold;
1278  Real better = betterThreshold(last);
1279  while (better < lastThreshold)
1280  {
1281  last = better;
1282  better = betterThreshold(last);
1283  }
1284  lastThreshold = last;
1285 
1286  // we reset the minimum stability (which might have been decreased below) to ensure that the increased sparsity
1287  // does not hurt the stability
1288  minStability = 2 * MINSTABILITY;
1289  }
1290 
1291  u.row.list.idx = thedim;
1292  u.row.start[thedim] = 0;
1293  u.row.max[thedim] = 0;
1294  u.row.len[thedim] = 0;
1295 
1296  u.col.list.idx = thedim;
1297  u.col.start[thedim] = 0;
1298  u.col.max[thedim] = 0;
1299  u.col.len[thedim] = 0;
1300 
1301  for (;;)
1302  {
1303  ///@todo if the factorization fails with stat = SINGULAR, distinuish between proven singularity (e.g., because of
1304  ///an empty column) and singularity due to numerics, that could be avoided by changing minStability and
1305  ///lastThreshold; in the first case, we want to abort, otherwise change the numerics
1306  stat = OK;
1307  factor(matrix, lastThreshold, epsilon);
1308 
1309  // finish if the factorization is stable
1310  if (stability() >= minStability)
1311  break;
1312 
1313  // otherwise, we increase the Markowitz threshold
1314  Real x = lastThreshold;
1316 
1317  // until it doesn't change anymore at its maximum value
1318  if (EQ(x, lastThreshold))
1319  break;
1320 
1321  // we relax the stability requirement
1322  minStability /= 2.0;
1323 
1324  MSG_INFO3( (*spxout), (*spxout) << "ISLUFA01 refactorizing with increased Markowitz threshold: "
1325  << lastThreshold << std::endl; )
1326  }
1327  MSG_DEBUG( std::cout << "DSLUFA02 threshold = " << lastThreshold
1328  << "\tstability = " << stability()
1329  << "\tminStability = " << minStability << std::endl; )
1330  MSG_DEBUG(
1331  int i;
1332  FILE* fl = fopen("dump.lp", "w");
1333  std::cout << "DSLUFA03 Basis:\n";
1334  int j = 0;
1335  for (i = 0; i < dim(); ++i)
1336  j += matrix[i]->size();
1337  for (i = 0; i < dim(); ++i)
1338  {
1339  for (j = 0; j < matrix[i]->size(); ++j)
1340  fprintf(fl, "%8d %8d %16g\n",
1341  i + 1, matrix[i]->index(j) + 1, matrix[i]->value(j));
1342  }
1343  fclose(fl);
1344  std::cout << "DSLUFA04 LU-Factors:" << std::endl;
1345  dump();
1346 
1347  std::cout << "DSLUFA05 threshold = " << lastThreshold
1348  << "\tstability = " << stability() << std::endl;
1349  )
1350 
1351  assert(isConsistent());
1352  return Status(stat);
1353 }
1354 
1355 
1357 {
1358 #ifdef ENABLE_CONSISTENCY_CHECKS
1359  return CLUFactor::isConsistent();
1360 #else
1361  return true;
1362 #endif
1363 }
1364 
1365 void SLUFactor::dump() const
1366 {
1367  CLUFactor::dump();
1368 }
1369 } // namespace soplex
int * len
used nonzeros per row vectors
Definition: clufactor.h:128
int updateType
type of updates to be used.
Definition: clufactor.h:167
const R * values() const
Returns array values.
Definition: ssvectorbase.h:288
DVector vec
Temporary vector.
Definition: slufactor.h:63
void updateNoClear(int p_col, const Real *p_work, const int *p_idx, int num)
Definition: clufactor.cpp:1320
Dring list
Double linked ringlist of vector indices in the order they appear in the column file.
Definition: clufactor.h:137
void solve2right4update(SSVector &x, Vector &y, const SVector &b, SSVector &d)
Solves and .
Definition: slufactor.cpp:106
Timer::TYPE timerType
Definition: slufactor.h:90
void reDim(int newdim, const bool setZero=true)
Resets DVectorBase&#39;s dimension to newdim.
Definition: dvectorbase.h:249
int * max
maximum available nonzeros per row: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]...
Definition: clufactor.h:129
SPxOut * spxout
message handler
Definition: slinsolver.h:209
#define REAL(x)
Definition: spxdefines.h:203
int vSolveRight4update(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *forest, int *forestNum, int *forestIdx)
Definition: clufactor.cpp:5637
Status load(const SVector *vec[], int dim)
Definition: slufactor.cpp:1224
Memory allocation routines.
void vSolveRight4update3sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real eps2, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2, Real eps3, Real *vec3, int *idx3, Real *rhs3, int *ridx3, int &rn3, Real *forest, int *forestNum, int *forestIdx)
sparse version of above method
Definition: clufactor.cpp:6055
SSVectorBase< R > & assign(const SVectorBase< S > &rhs)
Assigns only the elements of rhs.
Definition: basevectors.h:818
UpdateType uptype
the current UpdateType.
Definition: slufactor.h:73
SSVector ssvec
Temporary semi-sparse vector.
Definition: slufactor.h:64
void dump() const
Definition: clufactor.cpp:2829
struct soplex::CLUFactor::U::Col col
void solveLeft(Vector &x, const Vector &b)
sparse version of solving one system of equations with transposed basis matrix
Definition: slufactor.cpp:363
Real epsilon
|x| < epsililon is considered to be 0.
Definition: slufactor.h:87
int * start
starting positions in val and idx
Definition: clufactor.h:147
Real * rval
values of rows of L
Definition: clufactor.h:173
bool isConsistent() const
Definition: clufactor.cpp:2951
Exception classes for SoPlex.
#define initDR(ring)
Definition: cring.h:23
Real minStability
minimum stability to achieve by setting threshold.
Definition: slufactor.h:85
Real maxabs
maximum abs number in L and U
Definition: clufactor.h:189
bool LT(Real a, Real b, Real eps=Param::epsilon())
returns true iff a < b + eps
Definition: spxdefines.h:383
virtual ~Timer()
Definition: timer.h:124
int vSolveRight4update3(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real eps2, Real *rhs2, int *ridx2, int rn2, Real *vec3, Real eps3, Real *rhs3, int *ridx3, int rn3, Real *forest, int *forestNum, int *forestIdx)
Definition: clufactor.cpp:5905
void unSetup()
Makes SSVectorBase not setup.
Definition: ssvectorbase.h:126
int size() const
Number of used indices.
Definition: svectorbase.h:152
U u
U matrix.
Definition: clufactor.h:200
Implementation of Sparse Linear Solver.
virtual ~SLUFactor()
destructor.
Definition: slufactor.cpp:1198
void solveRight4update(SSVector &x, const SVector &b)
Solves .
Definition: slufactor.cpp:66
static Real epsilonFactorization()
Definition: spxdefines.h:262
void vSolveLeft3sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2, Real *vec3, int *idx3, Real *rhs3, int *ridx3, int &rn3)
sparse version of solving 3 systems of equations
Definition: clufactor.cpp:6324
Real lastThreshold
pivoting threshold of last factorization
Definition: slufactor.h:76
int size
size of array idx
Definition: clufactor.h:141
R & value(int n)
Reference to value of n &#39;th nonzero.
Definition: svectorbase.h:254
int * rorig
original row permutation
Definition: clufactor.h:176
void update(int p_col, Real *p_work, const int *p_idx, int num)
Definition: clufactor.cpp:1276
L l
L matrix.
Definition: clufactor.h:198
Dring * elem
Array of ring elements.
Definition: clufactor.h:122
void solveLeft(Real *vec, Real *rhs)
Definition: clufactor.cpp:4248
Real initMaxabs
maximum abs number in initail Matrix
Definition: clufactor.h:188
Timer * solveTime
Time spent in solves.
Definition: slufactor.h:89
int vSolveLeft2(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real *rhs2, int *ridx2, int rn2)
Definition: clufactor.cpp:6232
Real * diag
Array of pivot elements.
Definition: clufactor.h:199
Sparse Linear Solver virtual base class.Class SLinSolver provides a class for solving sparse linear s...
Definition: slinsolver.h:43
Perm row
row permutation matrices
Definition: clufactor.h:195
void solveRight(Vector &x, const Vector &b)
Solves .
Definition: slufactor.cpp:41
int * altIndexMem()
Returns array indices.
Definition: ssvectorbase.h:300
Dring * elem
Array of ring elements.
Definition: clufactor.h:140
Perm col
column permutation matrices
Definition: clufactor.h:196
int * orig
orig[p] original index from p
Definition: clufactor.h:109
Real minThreshold
minimum threshold to use.
Definition: slufactor.h:83
Real lMemMult
factor of minimum Memory * number of nonzeros
Definition: clufactor.h:193
Wrapper for different output streams and verbosity levels.
void changeEta(int idx, SSVector &eta)
Definition: slufactor.cpp:601
Exception class for out of memory exceptions.This class is derived from the SoPlex exception base cla...
Definition: exceptions.h:70
virtual void start()=0
start timer, resume accounting user, system and real time.
Real * work
Working array: must always be left as 0!
Definition: clufactor.h:202
void forestUpdate(int col, Real *work, int num, int *nonz)
Performs the Forrest-Tomlin update of the LU factorization.
Definition: clufactor.cpp:710
virtual Real stop()=0
stop timer, return accounted user time.
void spx_alloc(T &p, int n=1)
Allocate memory.
Definition: spxalloc.h:48
R * altValues()
Returns array values.
Definition: ssvectorbase.h:307
R * get_ptr()
Conversion to C-style pointer.
Definition: vectorbase.h:403
int * len
used nonzeros per column vector
Definition: clufactor.h:148
SLUFactor()
default constructor.
Definition: slufactor.cpp:975
double Real
SOPLEX_DEBUG.
Definition: spxdefines.h:200
int solveCount
Number of solves.
Definition: slufactor.h:92
#define MSG_DEBUG(x)
Definition: spxdefines.h:127
int getFactorCount() const
number of factorizations performed
Definition: slufactor.h:236
int getSolveCount() const
number of solves performed
Definition: slufactor.h:251
int size
size of arrays val and idx
Definition: clufactor.h:159
SLinSolver::Status stat
Status indicator.
Definition: clufactor.h:184
int used
used entries of arrays idx and val
Definition: clufactor.h:124
void solveLright(Real *vec)
Definition: clufactor.cpp:3319
int startSize
size of array start
Definition: clufactor.h:162
bool isConsistent() const
consistency check.
Definition: slufactor.cpp:1356
int factorCount
Number of factorizations.
Definition: clufactor.h:205
Real * val
hold nonzero values: this is only initialized in the end of the factorization with DEFAULT updates...
Definition: clufactor.h:144
void clear()
Clears vector.
Definition: ssvectorbase.h:572
int firstUnused
number of first unused L vector
Definition: clufactor.h:164
SLUFactor & operator=(const SLUFactor &old)
assignment operator.
Definition: slufactor.cpp:946
static Timer * createTimer(Timer::TYPE ttype)
create timers and allocate memory for them
Definition: timerfactory.h:43
void vSolveRight4update2sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real eps2, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2, Real *forest, int *forestNum, int *forestIdx)
sparse version of above method
Definition: clufactor.cpp:5827
int thedim
dimension of factorized matrix
Definition: clufactor.h:186
struct soplex::CLUFactor::U::Row row
int * start
starting positions in val and idx
Definition: clufactor.h:127
bool usetup
TRUE iff update vector has been setup.
Definition: slufactor.h:72
SSVector eta
Definition: slufactor.h:74
void reDim(int newdim)
Resets dimension to newdim.
Definition: ssvectorbase.h:538
Debugging, floating point type and parameter definitions.
void setSize(int n)
Sets number of nonzeros (thereby unSetup SSVectorBase).
Definition: ssvectorbase.h:553
int size
size of arrays val and idx
Definition: clufactor.h:123
Real getFactorTime() const
time spent in factorizations
Definition: slufactor.h:226
bool EQ(Real a, Real b, Real eps=Param::epsilon())
returns true iff |a-b| <= eps
Definition: spxdefines.h:371
Real getSolveTime() const
time spent in solves
Definition: slufactor.h:241
R * get_ptr()
Only used in slufactor.cpp.
Definition: ssvectorbase.h:99
void spx_realloc(T &p, int n)
Change amount of allocated memory.
Definition: spxalloc.h:79
int * idx
indices of L vectors
Definition: clufactor.h:161
void factor(const SVector **vec, Real threshold, Real eps)
epsilon for zero detection
Definition: clufactor.cpp:2772
int * perm
perm[i] permuted index from i
Definition: clufactor.h:110
Implementation of Sparse Linear Solver.This class implements a SLinSolver interface by using the spar...
Definition: slufactor.h:41
Implementation of sparse LU factorization.This class implements a sparse LU factorization with either...
Definition: clufactor.h:37
void setup_and_assign(SSVectorBase< S > &rhs)
Sets up rhs vector, and assigns it.
Definition: ssvectorbase.h:690
Everything should be within this namespace.
Status status() const
Definition: slufactor.h:171
void assign(const SLUFactor &old)
used to implement the assignment operator
Definition: slufactor.cpp:752
int vSolveLeft(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn)
Definition: clufactor.cpp:6204
void solve3right4update(SSVector &x, Vector &y, Vector &z, const SVector &b, SSVector &d, SSVector &e)
Solves , and .
Definition: slufactor.cpp:223
int firstUpdate
number of first update L vector
Definition: clufactor.h:163
int vSolveLeft3(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real *rhs2, int *ridx2, int rn2, Real *vec3, Real *rhs3, int *ridx3, int rn3)
Definition: clufactor.cpp:6288
Real rowMemMult
factor of minimum Memory * number of nonzeros
Definition: clufactor.h:191
Dring list
Double linked ringlist of vector indices in the order they appear in the row file.
Definition: clufactor.h:119
int vSolveRight4update2(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real eps2, Real *rhs2, int *ridx2, int rn2, Real *forest, int *forestNum, int *forestIdx)
Definition: clufactor.cpp:5704
Timer * factorTime
Time spent in factorizations.
Definition: clufactor.h:204
int dim() const
Definition: slufactor.h:156
void clear()
Set vector to 0.
Definition: vectorbase.h:219
void forceSetup()
Forces setup status.
Definition: ssvectorbase.h:162
int * idx
hold row indices of nonzeros
Definition: clufactor.h:143
std::string statistics() const
Definition: slufactor.cpp:590
int * ridx
indices of rows of L
Definition: clufactor.h:174
Real colMemMult
factor of minimum Memory * number of nonzeros
Definition: clufactor.h:192
#define MAXUPDATES
maximum nr. of factorization updates allowed before refactorization.
Definition: slufactor.h:33
int * max
maximum available nonzeros per colunn: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]...
Definition: clufactor.h:149
R getEpsilon() const
Returns the non-zero epsilon used.
Definition: ssvectorbase.h:104
#define MINSTABILITY
Definition: slufactor.cpp:39
int * row
column indices of L vectors
Definition: clufactor.h:166
void vSolveLeft2sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2)
sparse version of solving 2 systems of equations
Definition: clufactor.cpp:6261
const int * indexMem() const
Returns array indices.
Definition: ssvectorbase.h:282
int * rbeg
start of rows in rval and ridx
Definition: clufactor.h:175
void solveRight(Real *vec, Real *rhs)
Definition: clufactor.cpp:3581
int nzCnt
number of nonzeros in U
Definition: clufactor.h:187
Real * val
hold nonzero values
Definition: clufactor.h:125
int * rperm
original row permutation
Definition: clufactor.h:177
int * start
starting positions in val and idx
Definition: clufactor.h:165
VectorBase< R > & assign(const SVectorBase< S > &vec)
Assign values of vec.
Definition: basevectors.h:80
int * idx
hold column indices of nonzeros
Definition: clufactor.h:126
int used
used entries of array idx
Definition: clufactor.h:142
Status
status flags of the SLinSolver class.
Definition: slinsolver.h:51
#define MSG_INFO3(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO3.
Definition: spxdefines.h:117
No matrix has yet been loaded.
Definition: slinsolver.h:62
void dump() const
prints the LU factorization to stdout.
Definition: slufactor.cpp:1365
Status change(int idx, const SVector &subst, const SSVector *eta=0)
Definition: slufactor.cpp:610
SSVector forest
? Update vector set up by solveRight4update() and solve2right4update()
Definition: slufactor.h:75
Wrapper for the system time query methods.
Definition: timer.h:76
Real * val
values of L vectors
Definition: clufactor.h:160
void spx_free(T &p)
Release memory.
Definition: spxalloc.h:109
static Real betterThreshold(Real th)
Definition: slufactor.cpp:1208
int size() const
Returns the number of nonzeros.
Definition: ssvectorbase.h:199
SLinSolver::Status Status
for convenience
Definition: slufactor.h:55
Real stability() const
Definition: slufactor.cpp:578