clufactor.cpp

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1996-2019 Konrad-Zuse-Zentrum                            */
/*                            fuer Informationstechnik Berlin                */
/*                                                                           */
/*  SoPlex is distributed under the terms of the ZIB Academic Licence.       */
/*                                                                           */
/*  You should have received a copy of the ZIB Academic License              */
/*  along with SoPlex; see the file COPYING. If not email to soplex@zib.de.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#include <assert.h>

#include "soplex/spxdefines.h"
#include "soplex/clufactor.h"
#include "soplex/cring.h"
#include "soplex/spxalloc.h"
#include "soplex/exceptions.h"
#include "soplex/stablesum.h"

namespace soplex
{

/* This number is used to decide wether a value is zero
 * or was explicitly set to zero.
 */
#define SOPLEX_FACTOR_MARKER     1e-100

static const Real verySparseFactor = 0.001;
static const Real verySparseFactor4right = 0.2;
static const Real verySparseFactor4left  = 0.1;

/* generic heap management */
static void enQueueMax(int* heap, int* size, int elem)
{
   int i, j;

   j = (*size)++;

   while(j > 0)
   {
      i = (j - 1) / 2;

      if(elem > heap[i])
      {
         heap[j] = heap[i];
         j = i;
      }
      else
         break;
   }

   heap[j] = elem;

#ifdef SOPLEX_DEBUG

   // no NDEBUG define, since this block is really expensive
   for(i = 1; i < *size; ++i)
      for(j = 0; j < i; ++j)
         assert(heap[i] != heap[j]);

#endif  /* SOPLEX_DEBUG */
}

static int deQueueMax(int* heap, int* size)
{
   int e, elem;
   int i, j, s;
   int e1, e2;

   elem = *heap;
   e = heap[s = --(*size)];
   --s;

   for(j = 0, i = 1; i < s; i = 2 * j + 1)
   {
      e1 = heap[i];
      e2 = heap[i + 1];

      if(e1 > e2)
      {
         if(e < e1)
         {
            heap[j] = e1;
            j = i;
         }
         else
         {
            heap[j] = e;
            return elem;
         }
      }
      else
      {
         if(e < e2)
         {
            heap[j] = e2;
            j = i + 1;
         }
         else
         {
            heap[j] = e;
            return elem;
         }
      }
   }

   if(i < *size && e < heap[i])
   {
      heap[j] = heap[i];
      j = i;
   }

   heap[j] = e;

   return elem;
}

static void enQueueMin(int* heap, int* size, int elem)
{
   int i, j;

   j = (*size)++;

   while(j > 0)
   {
      i = (j - 1) / 2;

      if(elem < heap[i])
      {
         heap[j] = heap[i];
         j = i;
      }
      else
         break;
   }

   heap[j] = elem;

#ifdef SOPLEX_DEBUG

   // no NDEBUG define, since this block is really expensive
   for(i = 1; i < *size; ++i)
      for(j = 0; j < i; ++j)
         assert(heap[i] != heap[j]);

#endif  /* SOPLEX_DEBUG */
}

static int deQueueMin(int* heap, int* size)
{
   int e, elem;
   int i, j, s;
   int e1, e2;

   elem = *heap;
   e = heap[s = --(*size)];
   --s;

   for(j = 0, i = 1; i < s; i = 2 * j + 1)
   {
      e1 = heap[i];
      e2 = heap[i + 1];

      if(e1 < e2)
      {
         if(e > e1)
         {
            heap[j] = e1;
            j = i;
         }
         else
         {
            heap[j] = e;
            return elem;
         }
      }
      else
      {
         if(e > e2)
         {
            heap[j] = e2;
            j = i + 1;
         }
         else
         {
            heap[j] = e;
            return elem;
         }
      }
   }

   if(i < *size && e > heap[i])
   {
      heap[j] = heap[i];
      j = i;
   }

   heap[j] = e;

   return elem;
}

/************************************************************/
CLUFactor::Temp::Temp()
   : s_mark(0)
   , s_max(0)
   , s_cact(0)
   , stage(0)
   , pivot_col(0)
   , pivot_colNZ(0)
   , pivot_row(0)
   , pivot_rowNZ(0)
{}

void CLUFactor::Temp::init(int p_dim)
{
   spx_realloc(s_max,  p_dim);
   spx_realloc(s_cact, p_dim);
   spx_realloc(s_mark, p_dim);
   stage = 0;
}

void CLUFactor::Temp::clear()
{
   if(s_mark != 0)
      spx_free(s_mark);

   if(s_cact != 0)
      spx_free(s_cact);

   if(s_max != 0)
      spx_free(s_max);

   if(pivot_col != 0)
      spx_free(pivot_col);

   if(pivot_colNZ != 0)
      spx_free(pivot_colNZ);

   if(pivot_row != 0)
      spx_free(pivot_row);

   if(pivot_rowNZ != 0)
      spx_free(pivot_rowNZ);
}

CLUFactor::Temp::~Temp()
{
   clear();
}

/************************************************************/
void CLUFactor::initPerm()
{

   for(int i = 0; i < thedim; ++i)
      row.orig[i] = row.perm[i] = col.orig[i] = col.perm[i] = -1;
}

/*****************************************************************************/

void CLUFactor::setPivot(const int p_stage,
                         const int p_col,
                         const int p_row,
                         const Real val)
{
   assert(row.perm[p_row] < 0);
   assert(col.perm[p_col] < 0);

   row.orig[p_stage] = p_row;
   col.orig[p_stage] = p_col;
   row.perm[p_row]   = p_stage;
   col.perm[p_col]   = p_stage;
   diag[p_row]       = REAL(1.0) / val;

   if(spxAbs(val) < Param::epsilonPivot())
   {
#ifndef NDEBUG
      MSG_ERROR(std::cerr
                << "LU pivot element is almost zero (< "
                << Param::epsilonPivot()
                << ") - Basis is numerically singular"
                << std::endl;
               )
#endif
      stat = SLinSolver::SINGULAR;
   }

   if(spxAbs(diag[p_row]) > maxabs)
      maxabs = spxAbs(diag[p_row]);
}

/*****************************************************************************/
/*
 *      Perform garbage collection on row file
 */
void CLUFactor::packRows()
{
   int n, i, j, l_row;
   Dring* ring, *list;

   int* l_ridx = u.row.idx;
   Real* l_rval = u.row.val;
   int* l_rlen = u.row.len;
   int* l_rmax = u.row.max;
   int* l_rbeg = u.row.start;

   n = 0;
   list = &(u.row.list);

   for(ring = list->next; ring != list; ring = ring->next)
   {
      l_row = ring->idx;

      if(l_rbeg[l_row] != n)
      {
         do
         {
            l_row = ring->idx;
            i = l_rbeg[l_row];
            assert(l_rlen[l_row] <= l_rmax[l_row]);
            l_rbeg[l_row] = n;
            l_rmax[l_row] = l_rlen[l_row];
            j = i + l_rlen[l_row];

            for(; i < j; ++i, ++n)
            {
               assert(n <= i);
               l_ridx[n] = l_ridx[i];
               l_rval[n] = l_rval[i];
            }

            ring = ring->next;
         }
         while(ring != list);

         goto terminatePackRows;
      }

      n += l_rlen[l_row];

      l_rmax[l_row] = l_rlen[l_row];
   }

terminatePackRows:

   u.row.max[thedim] = 0;
   u.row.used = n;
}

/*****************************************************************************/
/*
 *      Perform garbage collection on column file
 */
void CLUFactor::forestPackColumns()
{
   int n, i, j, colno;
   Dring* ring, *list;

   Real* cval = u.col.val;
   int* cidx = u.col.idx;
   int* clen = u.col.len;
   int* cmax = u.col.max;
   int* cbeg = u.col.start;

   n = 0;
   list = &u.col.list;

   for(ring = list->next; ring != list; ring = ring->next)
   {
      colno = ring->idx;

      if(cbeg[colno] != n)
      {
         do
         {
            colno = ring->idx;
            i = cbeg[colno];
            cbeg[colno] = n;
            cmax[colno] = clen[colno];
            j = i + clen[colno];

            for(; i < j; ++i)
            {
               cval[n] = cval[i];
               cidx[n++] = cidx[i];
            }

            ring = ring->next;
         }
         while(ring != list);

         goto terminatePackColumns;
      }

      n += clen[colno];

      cmax[colno] = clen[colno];
   }

terminatePackColumns :

   u.col.used = n;
   u.col.max[thedim] = 0;
}

/*
 *      Make row of fac large enough to hold len nonzeros.
 */
void CLUFactor::remaxRow(int p_row, int len)
{
   assert(u.row.max[p_row] < len);

   if(u.row.elem[p_row].next == &(u.row.list))      /* last in row file */
   {
      int delta = len - u.row.max[p_row];

      if(delta > u.row.size - u.row.used)
      {
         packRows();
         delta = len - u.row.max[p_row];  // packRows() changes u.row.max[] !

         if(u.row.size < rowMemMult * u.row.used + len)
            minRowMem(2 * u.row.used + len);

         /* minRowMem(rowMemMult * u.row.used + len); */
      }

      assert(delta <= u.row.size - u.row.used

             && "ERROR: could not allocate memory for row file");

      u.row.used += delta;
      u.row.max[p_row] = len;
   }
   else                        /* row must be moved to end of row file */
   {
      int i, j, k;
      int* idx;
      Real* val;
      Dring* ring;

      if(len > u.row.size - u.row.used)
      {
         packRows();

         if(u.row.size < rowMemMult * u.row.used + len)
            minRowMem(2 * u.row.used + len);

         /* minRowMem(rowMemMult * u.row.used + len);*/
      }

      assert(len <= u.row.size - u.row.used

             && "ERROR: could not allocate memory for row file");

      j = u.row.used;
      i = u.row.start[p_row];
      k = u.row.len[p_row] + i;
      u.row.start[p_row] = j;
      u.row.used += len;

      u.row.max[u.row.elem[p_row].prev->idx] += u.row.max[p_row];
      u.row.max[p_row] = len;
      removeDR(u.row.elem[p_row]);
      ring = u.row.list.prev;
      init2DR(u.row.elem[p_row], *ring);

      idx = u.row.idx;
      val = u.row.val;

      for(; i < k; ++i, ++j)
      {
         val[j] = val[i];
         idx[j] = idx[i];
      }
   }

   assert(u.row.start[u.row.list.prev->idx] + u.row.max[u.row.list.prev->idx]

          == u.row.used);
}

/*************************************************************************/
/*
 *      Perform garbage collection on column file
 */
void CLUFactor::packColumns()
{
   int n, i, j, l_col;
   Dring* ring, *list;

   int* l_cidx = u.col.idx;
   int* l_clen = u.col.len;
   int* l_cmax = u.col.max;
   int* l_cbeg = u.col.start;

   n = 0;
   list = &(u.col.list);

   for(ring = list->next; ring != list; ring = ring->next)
   {
      l_col = ring->idx;

      if(l_cbeg[l_col] != n)
      {
         do
         {
            l_col = ring->idx;
            i = l_cbeg[l_col];
            l_cbeg[l_col] = n;
            l_cmax[l_col] = l_clen[l_col];
            j = i + l_clen[l_col];

            for(; i < j; ++i)
               l_cidx[n++] = l_cidx[i];

            ring = ring->next;
         }
         while(ring != list);

         goto terminatePackColumns;
      }

      n += l_clen[l_col];

      l_cmax[l_col] = l_clen[l_col];
   }

terminatePackColumns :

   u.col.used = n;
   u.col.max[thedim] = 0;
}

/*
 *      Make column col of fac large enough to hold len nonzeros.
 */
void CLUFactor::remaxCol(int p_col, int len)
{
   assert(u.col.max[p_col] < len);

   if(u.col.elem[p_col].next == &(u.col.list))      /* last in column file */
   {
      int delta = len - u.col.max[p_col];

      if(delta > u.col.size - u.col.used)
      {
         packColumns();
         delta = len - u.col.max[p_col];

         if(u.col.size < colMemMult * u.col.used + len)
            minColMem(2 * u.col.used + len);

         /* minColMem(colMemMult * u.col.used + len); */
      }

      assert(delta <= u.col.size - u.col.used

             && "ERROR: could not allocate memory for column file");

      u.col.used += delta;
      u.col.max[p_col] = len;
   }

   else                        /* column must be moved to end of column file */
   {
      int i, j, k;
      int* idx;
      Dring* ring;

      if(len > u.col.size - u.col.used)
      {
         packColumns();

         if(u.col.size < colMemMult * u.col.used + len)
            minColMem(2 * u.col.used + len);

         /* minColMem(colMemMult * u.col.used + len); */
      }

      assert(len <= u.col.size - u.col.used

             && "ERROR: could not allocate memory for column file");

      j = u.col.used;
      i = u.col.start[p_col];
      k = u.col.len[p_col] + i;
      u.col.start[p_col] = j;
      u.col.used += len;

      u.col.max[u.col.elem[p_col].prev->idx] += u.col.max[p_col];
      u.col.max[p_col] = len;
      removeDR(u.col.elem[p_col]);
      ring = u.col.list.prev;
      init2DR(u.col.elem[p_col], *ring);

      idx = u.col.idx;

      for(; i < k; ++i)
         idx[j++] = idx[i];
   }
}

/*
 *      Make column col of fac large enough to hold len nonzeros.
 */
void CLUFactor::forestReMaxCol(int p_col, int len)
{
   assert(u.col.max[p_col] < len);

   if(u.col.elem[p_col].next == &(u.col.list))      /* last in column file */
   {
      int delta = len - u.col.max[p_col];

      if(delta > u.col.size - u.col.used)
      {
         forestPackColumns();
         delta = len - u.col.max[p_col];

         if(u.col.size < colMemMult * u.col.used + len)
            forestMinColMem(int(colMemMult * u.col.used + len));
      }

      assert(delta <= u.col.size - u.col.used

             && "ERROR: could not allocate memory for column file");

      u.col.used += delta;
      u.col.max[p_col] = len;
   }

   else                        /* column must be moved to end of column file */
   {
      int i, j, k;
      int* idx;
      Real* val;
      Dring* ring;

      if(len > u.col.size - u.col.used)
      {
         forestPackColumns();

         if(u.col.size < colMemMult * u.col.used + len)
            forestMinColMem(int(colMemMult * u.col.used + len));
      }

      assert(len <= u.col.size - u.col.used

             && "ERROR: could not allocate memory for column file");

      j = u.col.used;
      i = u.col.start[p_col];
      k = u.col.len[p_col] + i;
      u.col.start[p_col] = j;
      u.col.used += len;

      u.col.max[u.col.elem[p_col].prev->idx] += u.col.max[p_col];
      u.col.max[p_col] = len;
      removeDR(u.col.elem[p_col]);
      ring = u.col.list.prev;
      init2DR(u.col.elem[p_col], *ring);

      idx = u.col.idx;
      val = u.col.val;

      for(; i < k; ++i)
      {
         val[j] = val[i];
         idx[j++] = idx[i];
      }
   }
}

/*****************************************************************************/

/**
 *   \brief Performs the Forrest-Tomlin update of the LU factorization.
 *
 *   BH: I suppose this is implemented as described in UH Suhl, LM Suhl: A fast LU
 *       update for linear programming, Annals of OR 43, p. 33-47, 1993.
 *
 *   @param  p_col      Index of basis column to replace.
 *   @param  p_work     Dense vector to substitute in the basis.
 *   @param  num        Number of nonzeros in vector represented by p_work.
 *   @param  nonz       Indices of nonzero elements in vector p_work.
 *
 *   The parameters num and nonz are used for the following optimization: If both
 *   are nonzero, indices of the nonzero elements provided in nonz (num giving
 *   their number) allow to access only those nonzero entries.  Otherwise we have
 *   to go through the entire dense vector element by element.
 *
 *   After copying p_work into U, p_work is used to expand the row r, which is
 *   needed to restore the triangular structure of U.
 *
 *   Also num and nonz are used to maintain a heap if there are only very few
 *   nonzeros to be eliminated. This is plainly wrong if the method is called with
 *   nonz==0, see todo at the corresponding place below.
 *
 *   @throw SPxStatusException if the loaded matrix is singular
 *
 *   @todo Use an extra member variable as a buffer for working with the dense
 *         row instead of misusing p_work. I think that should be as efficient and
 *         much cleaner.
 */

void CLUFactor::forestUpdate(int p_col, Real* p_work, int num, int* nonz)
{
   int i, j, k, h, m, n;
   int ll, c, r, rowno;
   Real x;

   Real* lval;
   int* lidx;
   int* lbeg = l.start;

   Real* cval;
   int* cidx = u.col.idx;
   int* cmax = u.col.max;
   int* clen = u.col.len;
   int* cbeg = u.col.start;

   Real* rval = u.row.val;
   int* ridx = u.row.idx;
   int* rmax = u.row.max;
   int* rlen = u.row.len;
   int* rbeg = u.row.start;

   int* rperm = row.perm;
   int* rorig = row.orig;
   int* cperm = col.perm;
   int* corig = col.orig;

   Real l_maxabs = maxabs;
   int dim = thedim;

   /*  Remove column p_col from U
    */
   j = cbeg[p_col];
   i = clen[p_col];
   nzCnt -= i;

   for(i += j - 1; i >= j; --i)
   {
      m = cidx[i];          // remove column p_col from row m
      k = rbeg[m];
      h = --(rlen[m]) + k;    // decrease length of row m

      while(ridx[k] != p_col)
         ++k;

      assert(k <= h);       // k is the position of p_col, h is last position

      ridx[k] = ridx[h];    // store last index at the position of p_col

      rval[k] = rval[h];
   }

   /*  Insert new vector column p_col thereby determining the highest permuted
    *       row index r.
    *
    *       Distinguish between optimized call (num > 0, nonz != 0) and
    *       non-optimized one.
    */
   assert(num);   // otherwise the assert( nonz != 0 ) below should fail

   if(num)
   {
      // Optimized call.
      assert(nonz != 0);

      clen[p_col] = 0;

      if(num > cmax[p_col])
         forestReMaxCol(p_col, num);

      cidx = u.col.idx;

      cval = u.col.val;

      k = cbeg[p_col];

      r = 0;

      for(j = 0; j < num; ++j)
      {
         i = nonz[j];
         x = p_work[i];
         p_work[i] = 0.0;

         if(isNotZero(x, Param::epsilonUpdate()))
         {
            if(spxAbs(x) > l_maxabs)
               l_maxabs = spxAbs(x);

            /* insert to column file */
            assert(k - cbeg[p_col] < cmax[p_col]);

            cval[k] = x;

            cidx[k++] = i;

            /* insert to row file */
            if(rmax[i] <= rlen[i])
            {
               remaxRow(i, rlen[i] + 1);
               rval = u.row.val;
               ridx = u.row.idx;
            }

            h = rbeg[i] + (rlen[i])++;

            rval[h] = x;
            ridx[h] = p_col;

            /* check permuted row index */

            if(rperm[i] > r)
               r = rperm[i];
         }
      }

      nzCnt += (clen[p_col] = k - cbeg[p_col]);
   }
   else
   {
      // Non-optimized call: We have to access all elements of p_work.
      assert(nonz == 0);

      /*
       *      clen[col] = 0;
       *      reMaxCol(fac, col, dim);
       */
      cidx = u.col.idx;
      cval = u.col.val;
      k = cbeg[p_col];
      j = k + cmax[p_col];
      r = 0;

      for(i = 0; i < dim; ++i)
      {
         x = p_work[i];
         p_work[i] = 0.0;

         if(isNotZero(x, Param::epsilonUpdate()))
         {
            if(spxAbs(x) > l_maxabs)
               l_maxabs = spxAbs(x);

            /* insert to column file */
            if(k >= j)
            {
               clen[p_col] = k - cbeg[p_col];
               forestReMaxCol(p_col, dim - i);
               cidx = u.col.idx;
               cval = u.col.val;
               k = cbeg[p_col];
               j = k + cmax[p_col];
               k += clen[p_col];
            }

            assert(k - cbeg[p_col] < cmax[p_col]);

            cval[k] = x;
            cidx[k++] = i;

            /* insert to row file */

            if(rmax[i] <= rlen[i])
            {
               remaxRow(i, rlen[i] + 1);
               rval = u.row.val;
               ridx = u.row.idx;
            }

            h = rbeg[i] + (rlen[i])++;

            rval[h] = x;
            ridx[h] = p_col;

            /* check permuted row index */

            if(rperm[i] > r)
               r = rperm[i];
         }
      }

      nzCnt += (clen[p_col] = k - cbeg[p_col]);

      if(cbeg[p_col] + cmax[p_col] == u.col.used)
      {
         u.col.used -= cmax[p_col];
         cmax[p_col] = clen[p_col];
         u.col.used += cmax[p_col];
      }
   }

   c = cperm[p_col];

   if(r > c)                          /* Forest Tomlin update */
   {
      /*      update permutations
       */
      j = rorig[c];

      // memmove is more efficient than a for loop
      // for ( i = c; i < r; ++i )
      //    rorig[i] = rorig[i + 1];
      memmove(&rorig[c], &rorig[c + 1], (unsigned int)(r - c) * sizeof(int));

      rorig[r] = j;

      for(i = c; i <= r; ++i)
         rperm[rorig[i]] = i;

      j = corig[c];

      // memmove is more efficient than a for loop
      // for ( i = c; i < r; ++i )
      //    corig[i] = corig[i + 1];
      memmove(&corig[c], &corig[c + 1], (unsigned int)(r - c) * sizeof(int));

      corig[r] = j;

      for(i = c; i <= r; ++i)
         cperm[corig[i]] = i;


      rowno = rorig[r];

      j = rbeg[rowno];

      i = rlen[rowno];

      nzCnt -= i;

      if(i < verySparseFactor * (dim - c))      // few nonzeros to be eliminated
      {
         /**
          *          The following assert is obviously violated if this method is called
          *          with nonzero==0.
          *
          *          @todo Use an extra member variable as a buffer for the heap instead of
          *                misusing nonz and num. The method enQueueMin() seems to
          *                sort the nonzeros or something, for which it only needs
          *                some empty vector of size num.
          */
         assert(nonz != 0);

         /*  move row r from U to p_work
          */
         num = 0;

         for(i += j - 1; i >= j; --i)
         {
            k = ridx[i];
            p_work[k] = rval[i];
            enQueueMin(nonz, &num, cperm[k]);
            m = --(clen[k]) + cbeg[k];

            for(h = m; cidx[h] != rowno; --h)
               ;

            cidx[h] = cidx[m];

            cval[h] = cval[m];
         }


         /*  Eliminate row r from U to L file
          */
         ll = makeLvec(r - c, rowno);

         lval = l.val;

         lidx = l.idx;

         assert((num == 0) || (nonz != 0));

         /* for(i = c; i < r; ++i)       */
         while(num)
         {
#ifndef NDEBUG
            // The numbers seem to be often 1e-100, is this ok ?

            for(i = 0; i < num; ++i)
               assert(p_work[corig[nonz[i]]] != 0.0);

#endif // NDEBUG
            i = deQueueMin(nonz, &num);

            if(i == r)
               break;

            k = corig[i];

            assert(p_work[k] != 0.0);

            n = rorig[i];

            x = p_work[k] * diag[n];

            lidx[ll] = n;

            lval[ll] = x;

            p_work[k] = 0.0;

            ll++;

            if(spxAbs(x) > l_maxabs)
               l_maxabs = spxAbs(x);

            j = rbeg[n];

            m = rlen[n] + j;

            for(; j < m; ++j)
            {
               int jj = ridx[j];
               Real y = p_work[jj];

               if(y == 0)
                  enQueueMin(nonz, &num, cperm[jj]);

               y -= x * rval[j];

               p_work[jj] = y + ((y == 0) ? SOPLEX_FACTOR_MARKER : 0.0);
            }
         }

         if(lbeg[l.firstUnused - 1] == ll)
            (l.firstUnused)--;
         else
            lbeg[l.firstUnused] = ll;


         /*  Set diagonal value
          */
         if(i != r)
         {
            stat = SLinSolver::SINGULAR;
            throw SPxStatusException("XFORE01 The loaded matrix is singular");
         }

         k = corig[r];

         x = p_work[k];
         diag[rowno] = 1 / x;
         p_work[k] = 0.0;


         /*  make row large enough to fit all nonzeros.
          */

         if(rmax[rowno] < num)
         {
            rlen[rowno] = 0;
            remaxRow(rowno, num);
            rval = u.row.val;
            ridx = u.row.idx;
         }

         nzCnt += num;

         /*  Insert work to updated row thereby clearing work;
          */
         n = rbeg[rowno];

         for(i = 0; i < num; ++i)
         {
            j = corig[nonz[i]];
            x = p_work[j];

            // BH 2005-08-24: This if is very important. It may well happen that
            // during the elimination of row r a nonzero elements cancels out
            // and becomes zero. This would lead to an infinite loop in the
            // above elimination code, since the corresponding column index would
            // be enqueued for further elimination again and agian.

            if(x != 0.0)
            {
               if(spxAbs(x) > l_maxabs)
                  l_maxabs = spxAbs(x);

               ridx[n] = j;

               rval[n] = x;

               p_work[j] = 0.0;

               ++n;

               if(clen[j] >= cmax[j])
               {
                  forestReMaxCol(j, clen[j] + 1);
                  cidx = u.col.idx;
                  cval = u.col.val;
               }

               cval[cbeg[j] + clen[j]] = x;

               cidx[cbeg[j] + clen[j]++] = rowno;
            }
         }

         rlen[rowno] = n - rbeg[rowno];
      }
      else            /* few nonzeros to be eliminated        */
      {
         /*  move row r from U to p_work
          */
         for(i += j - 1; i >= j; --i)
         {
            k = ridx[i];
            p_work[k] = rval[i];
            m = --(clen[k]) + cbeg[k];

            for(h = m; cidx[h] != rowno; --h)
               ;

            cidx[h] = cidx[m];

            cval[h] = cval[m];
         }


         /*  Eliminate row r from U to L file
          */
         ll = makeLvec(r - c, rowno);

         lval = l.val;

         lidx = l.idx;

         for(i = c; i < r; ++i)
         {
            k = corig[i];

            if(p_work[k] != 0.0)
            {
               n = rorig[i];
               x = p_work[k] * diag[n];
               lidx[ll] = n;
               lval[ll] = x;
               p_work[k] = 0.0;
               ll++;

               if(spxAbs(x) > l_maxabs)
                  l_maxabs = spxAbs(x);

               j = rbeg[n];

               m = rlen[n] + j;

               for(; j < m; ++j)
                  p_work[ridx[j]] -= x * rval[j];
            }
         }

         if(lbeg[l.firstUnused - 1] == ll)
            (l.firstUnused)--;
         else
            lbeg[l.firstUnused] = ll;


         /*  Set diagonal value
          */
         k = corig[r];

         x = p_work[k];

         if(x == 0.0)
         {
            stat = SLinSolver::SINGULAR;
            throw SPxStatusException("XFORE02 The loaded matrix is singular");
            //            return;
         }

         diag[rowno] = 1 / x;

         p_work[k] = 0.0;


         /*  count remaining nonzeros in work and make row large enough
          *  to fit them all.
          */
         n = 0;

         for(i = r + 1; i < dim; ++i)
            if(p_work[corig[i]] != 0.0)
               n++;

         if(rmax[rowno] < n)
         {
            rlen[rowno] = 0;
            remaxRow(rowno, n);
            rval = u.row.val;
            ridx = u.row.idx;
         }

         nzCnt += n;

         /*  Insert p_work to updated row thereby clearing p_work;
          */
         n = rbeg[rowno];

         for(i = r + 1; i < dim; ++i)
         {
            j = corig[i];
            x = p_work[j];

            if(x != 0.0)
            {
               if(spxAbs(x) > l_maxabs)
                  l_maxabs = spxAbs(x);

               ridx[n] = j;

               rval[n] = x;

               p_work[j] = 0.0;

               ++n;

               if(clen[j] >= cmax[j])
               {
                  forestReMaxCol(j, clen[j] + 1);
                  cidx = u.col.idx;
                  cval = u.col.val;
               }

               cval[cbeg[j] + clen[j]] = x;

               cidx[cbeg[j] + clen[j]++] = rowno;
            }
         }

         rlen[rowno] = n - rbeg[rowno];
      }
   }

   else if(r == c)
   {
      /*  Move diagonal element to diag.  Note, that it must be the last
       *  element, since it has just been inserted above.
       */
      rowno = rorig[r];
      i = rbeg[rowno] + --(rlen[rowno]);
      diag[rowno] = 1 / rval[i];

      for(j = i = --(clen[p_col]) + cbeg[p_col]; cidx[i] != rowno; --i)
         ;

      cidx[i] = cidx[j];

      cval[i] = cval[j];
   }
   else /* r < c */
   {
      stat = SLinSolver::SINGULAR;
      throw SPxStatusException("XFORE03 The loaded matrix is singular");
      //      return;
   }

   maxabs = l_maxabs;

   assert(isConsistent());
   stat = SLinSolver::OK;
}

void CLUFactor::update(int p_col, Real* p_work, const int* p_idx, int num)
{
   int ll, i, j;
   int* lidx;
   Real* lval;
   Real x, rezi;

   assert(p_work[p_col] != 0.0);
   rezi = 1 / p_work[p_col];
   p_work[p_col] = 0.0;

   ll = makeLvec(num, p_col);
   //   ll = fac->makeLvec(num, col);
   lval = l.val;
   lidx = l.idx;

   for(i = num - 1; (j = p_idx[i]) != p_col; --i)
   {
      lidx[ll] = j;
      lval[ll] = rezi * p_work[j];
      p_work[j] = 0.0;
      ++ll;
   }

   lidx[ll] = p_col;

   lval[ll] = 1 - rezi;
   ++ll;

   for(--i; i >= 0; --i)
   {
      j = p_idx[i];
      lidx[ll] = j;
      lval[ll] = x = rezi * p_work[j];
      p_work[j] = 0.0;
      ++ll;

      if(spxAbs(x) > maxabs)
         maxabs = spxAbs(x);
   }

   stat = SLinSolver::OK;
}

void CLUFactor::updateNoClear(
   int p_col,
   const Real* p_work,
   const int* p_idx,
   int num)
{
   int ll, i, j;
   int* lidx;
   Real* lval;
   Real x, rezi;

   assert(p_work[p_col] != 0.0);
   rezi = 1 / p_work[p_col];
   ll = makeLvec(num, p_col);
   //ll = fac->makeLvec(num, col);
   lval = l.val;
   lidx = l.idx;

   for(i = num - 1; (j = p_idx[i]) != p_col; --i)
   {
      lidx[ll] = j;
      lval[ll] = rezi * p_work[j];
      ++ll;
   }

   lidx[ll] = p_col;

   lval[ll] = 1 - rezi;
   ++ll;

   for(--i; i >= 0; --i)
   {
      j = p_idx[i];
      lidx[ll] = j;
      lval[ll] = x = rezi * p_work[j];
      ++ll;

      if(spxAbs(x) > maxabs)
         maxabs = spxAbs(x);
   }

   stat = SLinSolver::OK;
}

/*****************************************************************************/
/*
 *      Temporary data structures.
 */

/*
 *        For the i=th column the situation might look like this:
 *
 * \begin{verbatim}
 *        idx     = ....................iiiIIIIII-----..............
 *        cbeg[i] =                     ^
 *        cact[i] =                        +----+
 *        clen[i] =                     +-------+
 *        cmax[i] =                     +------------+
 *
 *        Indices clen[i]-cbeg[i]:      ^^^
 * \end{verbatim}
 *        belong to column i, but have already been pivotal and don't belong to
 *        the active submatrix.
 */

/****************************************************************************/
/*
 *      Initialize row and column file of working matrix and
 *      mark column singletons.
 */
void CLUFactor::initFactorMatrix(const SVector** vec, const Real eps)
{

   Real x;
   int m;
   int tot;
   Dring* rring, *lastrring;
   Dring* cring, *lastcring;
   const SVector* psv;
   int* sing = temp.s_mark;

   /*  Initialize:
    *  - column file thereby remembering column singletons in |sing|.
    *  - nonzeros counts per row
    *  - total number of nonzeros
    */

   for(int i = 0; i < thedim; i++)
      u.row.max[i] = u.row.len[i] = 0;

   tot = 0;

   for(int i = 0; i < thedim; i++)
   {
      int k;

      psv = vec[i];
      k = psv->size();

      if(k > 1)
      {
         tot += k;

         for(int j = 0; j < k; ++j)
            u.row.max[psv->index(j)]++;
      }
      else if(k == 0)
      {
         stat = SLinSolver::SINGULAR;
         return;
      }
   }

   /*  Resize nonzero memory if necessary
    */
   minRowMem(int(rowMemMult * tot));

   minColMem(int(colMemMult * tot));

   minLMem(int(lMemMult * tot));


   /*  Initialize:
    *  - row ring lists
    *  - row vectors in file
    *  - column ring lists
    */
   u.row.start[0] = 0;

   rring = u.row.elem;

   lastrring = &(u.row.list);

   lastrring->idx = thedim;

   lastrring->next = rring;

   cring = u.col.elem;

   lastcring = &(u.col.list);

   lastcring->idx = thedim;

   lastcring->next = cring;

   m = 0;

   for(int i = 0; i < thedim; i++)
   {
      u.row.start[i] = m;
      m += u.row.max[i];

      rring->idx = i;
      rring->prev = lastrring;
      lastrring->next = rring;
      lastrring = rring;
      ++rring;

      cring->idx = i;
      cring->prev = lastcring;
      lastcring->next = cring;
      lastcring = cring;
      ++cring;
   }

   u.row.start[thedim]       = 0;

   u.row.max[thedim]       = 0;
   u.row.used = m;

   lastrring->next = &(u.row.list);
   lastrring->next->prev = lastrring;

   lastcring->next = &(u.col.list);
   lastcring->next->prev = lastcring;

   /*  Copy matrix to row and column file
    *  excluding and marking column singletons!
    */
   m = 0;
   temp.stage = 0;

   initMaxabs = 0;

   for(int i = 0; i < thedim; i++)
   {
      int nnonzeros;

      psv = vec[i];
      u.col.start[i] = m;

      /* check whether number of nonzeros above tolerance is 0, 1 or >= 2 */
      nnonzeros = 0;

      for(int j = 0; j < psv->size() && nnonzeros <= 1; j++)
      {
         if(isNotZero(psv->value(j), eps))
            nnonzeros++;
      }

      /* basis is singular due to empty column */
      if(nnonzeros == 0)
      {
         stat = SLinSolver::SINGULAR;
         return;
      }

      /* exclude column singletons */
      else if(nnonzeros == 1)
      {
         int j;

         /* find nonzero */

         for(j = 0; isZero(psv->value(j), eps); j++)
            ;

         assert(j < psv->size());

         /* basis is singular due to two linearly dependent column singletons */
         if(row.perm[psv->index(j)] >= 0)
         {
            stat = SLinSolver::SINGULAR;
            return;
         }

         /* update maximum absolute nonzero value */
         x = psv->value(j);

         if(spxAbs(x) > initMaxabs)
            initMaxabs = spxAbs(x);

         /* permute to front and mark as singleton */
         setPivot(temp.stage, i, psv->index(j), x);

         sing[temp.stage] = i;

         temp.stage++;

         /* set column length to zero */
         temp.s_cact[i] = u.col.len[i] = u.col.max[i] = 0;
      }

      /* add to active matrix if not a column singleton */
      else
      {
         int end;
         int k;

         /* go through all nonzeros in column */
         assert(nnonzeros >= 2);
         nnonzeros = 0;

         for(int j = 0; j < psv->size(); j++)
         {
            x = psv->value(j);

            if(isNotZero(x, eps))
            {
               /* add to column array */
               k = psv->index(j);
               u.col.idx[m] = k;
               m++;

               /* add to row array */
               end = u.row.start[k] + u.row.len[k];
               u.row.idx[end] = i;
               u.row.val[end] = x;
               u.row.len[k]++;

               /* update maximum absolute nonzero value */

               if(spxAbs(x) > initMaxabs)
                  initMaxabs = spxAbs(x);

               nnonzeros++;
            }
         }

         assert(nnonzeros >= 2);

         /* set column length */
         temp.s_cact[i] = u.col.len[i] = u.col.max[i] = nnonzeros;
      }
   }

   u.col.used = m;
}



/*****************************************************************************/
/*
 *      Remove column singletons
 */

void CLUFactor::colSingletons()
{
   int i, j, k, n;
   int len;
   int p_col, p_row, newrow;
   int* idx;
   int* rorig = row.orig;
   int* rperm = row.perm;
   int* sing = temp.s_mark;


   /*  Iteratively update column counts due to removed column singletons
    *  thereby removing new arising columns singletons
    *  and computing the index of the first row singleton (-1)
    *  until no more can be found.
    */

   for(i = 0; i < temp.stage; ++i)
   {
      p_row = rorig[i];
      assert(p_row >= 0);
      idx = &(u.row.idx[u.row.start[p_row]]);
      len = u.row.len[p_row];

      for(j = 0; j < len; ++j)
      {
         /*  Move pivotal nonzeros to front of column.
          */
         p_col = idx[j];
         assert(temp.s_cact[p_col] > 0);

         n = u.col.start[p_col] + u.col.len[p_col] - temp.s_cact[p_col];

         for(k = n; u.col.idx[k] != p_row; ++k)
            ;

         assert(k < u.col.start[p_col] + u.col.len[p_col]);

         u.col.idx[k] = u.col.idx[n];

         u.col.idx[n] = p_row;

         n = --(temp.s_cact[p_col]);          /* column nonzeros of ACTIVE matrix */

         if(n == 1)                   /* Here is another singleton */
         {
            newrow = u.col.idx[--u.col.len[p_col] + u.col.start[p_col]];

            /*      Ensure, matrix not singular
             */

            if(rperm[newrow] >= 0)
            {
               stat = SLinSolver::SINGULAR;
               return;
            }

            /*      Find singleton in row.
             */
            n = u.row.start[newrow] + (--(u.row.len[newrow]));

            for(k = n; u.row.idx[k] != p_col; --k)
               ;

            /*      Remove singleton from column.
             */
            setPivot(temp.stage, p_col, newrow, u.row.val[k]);

            sing[temp.stage++] = p_col;

            /*      Move pivot element to diag.
             */
            u.row.val[k] = u.row.val[n];

            u.row.idx[k] = u.row.idx[n];
         }
         else if(n == 0)
         {
            stat = SLinSolver::SINGULAR;
            return;
         }
      }
   }

   assert(temp.stage <= thedim);
}


/*****************************************************************************/
/*
 *      Remove row singletons
 */
void CLUFactor::rowSingletons()
{
   Real pval;
   int i, j, k, ll, r;
   int p_row, p_col, len, rs, lk;
   int* idx;
   int* rperm = row.perm;
   int* sing = temp.s_mark;

   /*  Mark row singletons
    */
   rs = temp.stage;

   for(i = 0; i < thedim; ++i)
   {
      if(rperm[i] < 0 && u.row.len[i] == 1)
         sing[temp.stage++] = i;
   }

   /*  Eliminate row singletons
    *  thereby marking newly arising ones
    *  until no more can be found.
    */
   for(; rs < temp.stage; ++rs)
   {
      /*      Move pivot element from row file to diag
       */
      p_row = sing[rs];
      j = u.row.start[p_row];
      p_col = u.row.idx[j];
      pval = u.row.val[j];
      setPivot(rs, p_col, p_row, pval);
      u.row.len[p_row] = 0;

      /*      Remove pivot column form workingmatrix
       *      thereby building up L vector.
       */
      idx = &(u.col.idx[u.col.start[p_col]]);
      i = temp.s_cact[p_col];                /* nr. nonzeros of new L vector */
      lk = makeLvec(i - 1, p_row);
      len = u.col.len[p_col];
      i = (u.col.len[p_col] -= i);         /* remove pivot column from U */

      for(; i < len; ++i)
      {
         r = idx[i];

         if(r != p_row)
         {
            /*      Find pivot column in row.
             */
            ll = --(u.row.len[r]);
            k = u.row.start[r] + ll;

            for(j = k; u.row.idx[j] != p_col; --j)
               ;

            assert(k >= u.row.start[r]);

            /*      Initialize L vector
             */
            l.idx[lk] = r;

            l.val[lk] = u.row.val[j] / pval;

            ++lk;

            /*      Remove pivot column from row.
             */
            u.row.idx[j] = u.row.idx[k];

            u.row.val[j] = u.row.val[k];

            /*      Check new row length.
             */
            if(ll == 1)
               sing[temp.stage++] = r;
            else if(ll == 0)
            {
               stat = SLinSolver::SINGULAR;
               return;
            }
         }
      }
   }
}


/*****************************************************************************/
/*
 *      Init nonzero number Ring lists
 *      and required entries of arrays max and mark
 */

void CLUFactor::initFactorRings()
{
   int i;
   int* rperm = row.perm;
   int* cperm = col.perm;
   CLUFactor::Pring* ring;

   assert(thedim >= 0);
   spx_alloc(temp.pivot_col,   thedim + 1);
   spx_alloc(temp.pivot_colNZ, thedim + 1);
   spx_alloc(temp.pivot_row,   thedim + 1);
   spx_alloc(temp.pivot_rowNZ, thedim + 1);

   for(i = thedim - temp.stage; i >= 0; --i)
   {
      initDR(temp.pivot_colNZ[i]);
      initDR(temp.pivot_rowNZ[i]);
   }

   for(i = 0; i < thedim; ++i)
   {
      if(rperm[i] < 0)
      {
         if(u.row.len[i] <= 0)
         {
            stat = SLinSolver::SINGULAR;
            return;
         }

         ring = &(temp.pivot_rowNZ[u.row.len[i]]);

         init2DR(temp.pivot_row[i], *ring);
         temp.pivot_row[i].idx = i;
         temp.s_max[i] = -1;
      }

      if(cperm[i] < 0)
      {
         if(temp.s_cact[i] <= 0)
         {
            stat = SLinSolver::SINGULAR;
            return;
         }

         ring = &(temp.pivot_colNZ[temp.s_cact[i]]);

         init2DR(temp.pivot_col[i], *ring);
         temp.pivot_col[i].idx = i;
         temp.s_mark[i] = 0;
      }
   }
}

void CLUFactor::freeFactorRings(void)
{

   if(temp.pivot_col)
      spx_free(temp.pivot_col);

   if(temp.pivot_colNZ)
      spx_free(temp.pivot_colNZ);

   if(temp.pivot_row)
      spx_free(temp.pivot_row);

   if(temp.pivot_rowNZ)
      spx_free(temp.pivot_rowNZ);
}


/*
 *      Eliminate all row singletons from nucleus.
 *      A row singleton may well be column singleton at the same time!
 */
void CLUFactor::eliminateRowSingletons()
{
   int i, j, k, ll, r;
   int len, lk;
   int pcol, prow;
   Real pval;
   int* idx;
   CLUFactor::Pring* sing;

   for(sing = temp.pivot_rowNZ[1].prev; sing != &(temp.pivot_rowNZ[1]); sing = sing->prev)
   {
      prow = sing->idx;
      i = u.row.start[prow];
      pcol = u.row.idx[i];
      pval = u.row.val[i];
      setPivot(temp.stage++, pcol, prow, pval);
      u.row.len[prow] = 0;
      removeDR(temp.pivot_col[pcol]);

      /*      Eliminate pivot column and build L vector.
       */
      i = temp.s_cact[pcol];

      if(i > 1)
      {
         idx = &(u.col.idx[u.col.start[pcol]]);
         len = u.col.len[pcol];
         lk = makeLvec(i - 1, prow);
         i = u.col.len[pcol] -= i;

         for(; (r = idx[i]) != prow; ++i)
         {
            /*      Find pivot column in row.
             */
            ll = --(u.row.len[r]);
            k = u.row.start[r] + ll;

            for(j = k; u.row.idx[j] != pcol; --j)
               ;

            assert(j >= u.row.start[r]);

            /*      Initialize L vector
             */
            l.idx[lk] = r;

            l.val[lk] = u.row.val[j] / pval;

            ++lk;

            /*      Remove pivot column from row.
             */
            u.row.idx[j] = u.row.idx[k];

            u.row.val[j] = u.row.val[k];

            /*      Move column to appropriate nonzero ring.
             */
            removeDR(temp.pivot_row[r]);

            init2DR(temp.pivot_row[r], temp.pivot_rowNZ[ll]);

            assert(row.perm[r] < 0);

            temp.s_max[r] = -1;
         }

         /* skip pivot element */
         assert(i < len && "ERROR: pivot column does not contain pivot row");

         for(++i; i < len; ++i)
         {
            /*      Find pivot column in row.
             */
            r = idx[i];
            ll = --(u.row.len[r]);
            k = u.row.start[r] + ll;

            for(j = k; u.row.idx[j] != pcol; --j)
               ;

            assert(j >= u.row.start[r]);

            /*      Initialize L vector
             */
            l.idx[lk] = r;

            l.val[lk] = u.row.val[j] / pval;

            ++lk;

            /*      Remove pivot column from row.
             */
            u.row.idx[j] = u.row.idx[k];

            u.row.val[j] = u.row.val[k];

            /*      Move column to appropriate nonzero ring.
             */
            removeDR(temp.pivot_row[r]);

            init2DR(temp.pivot_row[r], temp.pivot_rowNZ[ll]);

            assert(row.perm[r] < 0);

            temp.s_max[r] = -1;
         }
      }
      else
         u.col.len[pcol] -= i;
   }

   initDR(temp.pivot_rowNZ[1]);    /* Remove all row singletons from list */
}



/*
 *      Eliminate all column singletons from nucleus.
 *      A column singleton must not be row singleton at the same time!
 */
void CLUFactor::eliminateColSingletons()
{
   int i, j, k, m, c;
   int pcol, prow;
   CLUFactor::Pring* sing;

   for(sing = temp.pivot_colNZ[1].prev;
         sing != &(temp.pivot_colNZ[1]);
         sing = sing->prev)
   {
      /*      Find pivot value
       */
      pcol = sing->idx;
      j = --(u.col.len[pcol]) + u.col.start[pcol];   /* remove pivot column */
      prow = u.col.idx[j];
      removeDR(temp.pivot_row[prow]);

      j = --(u.row.len[prow]) + u.row.start[prow];

      for(i = j; (c = u.row.idx[i]) != pcol; --i)
      {
         m = u.col.len[c] + u.col.start[c] - (temp.s_cact[c])--;

         for(k = m; u.col.idx[k] != prow; ++k)
            ;

         u.col.idx[k] = u.col.idx[m];

         u.col.idx[m] = prow;

         m = temp.s_cact[c];

         removeDR(temp.pivot_col[c]);

         init2DR(temp.pivot_col[c], temp.pivot_colNZ[m]);

         assert(col.perm[c] < 0);
      }

      /*      remove pivot element from pivot row
       */
      setPivot(temp.stage++, pcol, prow, u.row.val[i]);

      u.row.idx[i] = u.row.idx[j];

      u.row.val[i] = u.row.val[j];

      j = u.row.start[prow];

      for(--i; i >= j; --i)
      {
         c = u.row.idx[i];
         m = u.col.len[c] + u.col.start[c] - (temp.s_cact[c])--;

         for(k = m; u.col.idx[k] != prow; ++k)
            ;

         u.col.idx[k] = u.col.idx[m];

         u.col.idx[m] = prow;

         m = temp.s_cact[c];

         removeDR(temp.pivot_col[c]);

         init2DR(temp.pivot_col[c], temp.pivot_colNZ[m]);

         assert(col.perm[c] < 0);
      }
   }

   initDR(temp.pivot_colNZ[1]);    /* Remove all column singletons from list */
}

/*
 * No singletons available: Select pivot elements.
 */
void CLUFactor::selectPivots(Real threshold)
{
   int ii;
   int i;
   int j;
   int k;
   int ll = -1; // This value should never be used.
   int kk;
   int m;
   int count;
   int num;
   int rw = -1; // This value should never be used.
   int cl = -1; // This value should never be used.
   int len;
   int beg;
   Real l_maxabs;
   Real x = REAL(0.0);   // This value should never be used.
   int mkwtz;
   int candidates;

   candidates = thedim - temp.stage - 1;

   if(candidates > 4)
      candidates = 4;

   num = 0;

   count = 2;

   for(;;)
   {
      ii = -1;

      if(temp.pivot_rowNZ[count].next != &(temp.pivot_rowNZ[count]))
      {
         rw = temp.pivot_rowNZ[count].next->idx;
         beg = u.row.start[rw];
         len = u.row.len[rw] + beg - 1;

         /*  set l_maxabs to maximum absolute value in row
          *  (compute it if necessary).
          */

         if((l_maxabs = temp.s_max[rw]) < 0)
         {
            l_maxabs = spxAbs(u.row.val[len]);

            for(i = len - 1; i >= beg; --i)
               if(l_maxabs < spxAbs(u.row.val[i]))
                  l_maxabs = spxAbs(u.row.val[i]);

            temp.s_max[rw] = l_maxabs;               /* ##### */
         }

         l_maxabs *= threshold;

         /*  select pivot element with lowest markowitz number in row
          */
         mkwtz = thedim + 1;

         for(i = len; i >= beg; --i)
         {
            k = u.row.idx[i];
            j = temp.s_cact[k];
            x = u.row.val[i];

            if(j < mkwtz && spxAbs(x) > l_maxabs)
            {
               mkwtz = j;
               cl = k;
               ii = i;

               if(j <= count)               /* ##### */
                  break;
            }
         }
      }
      else if(temp.pivot_colNZ[count].next != &(temp.pivot_colNZ[count]))
      {
         cl = temp.pivot_colNZ[count].next->idx;
         beg = u.col.start[cl];
         len = u.col.len[cl] + beg - 1;
         beg = len - temp.s_cact[cl] + 1;
         assert(count == temp.s_cact[cl]);

         /*  select pivot element with lowest markowitz number in column
          */
         mkwtz = thedim + 1;

         for(i = len; i >= beg; --i)
         {
            k = u.col.idx[i];
            j = u.row.len[k];

            if(j < mkwtz)
            {
               /*  ensure that element (cl,k) is stable.
                */
               if(temp.s_max[k] > 0)
               {
                  /*  case 1: l_maxabs is known
                   */
                  for(m = u.row.start[k], kk = m + u.row.len[k] - 1;
                        kk >= m; --kk)
                  {
                     if(u.row.idx[kk] == cl)
                     {
                        x = u.row.val[kk];
                        ll = kk;
                        break;
                     }
                  }

                  l_maxabs = temp.s_max[k];
               }
               else
               {
                  /*  case 2: l_maxabs needs to be computed
                   */
                  m = u.row.start[k];
                  l_maxabs = spxAbs(u.row.val[m]);

                  for(kk = m + u.row.len[k] - 1; kk >= m; --kk)
                  {
                     if(l_maxabs < spxAbs(u.row.val[kk]))
                        l_maxabs = spxAbs(u.row.val[kk]);

                     if(u.row.idx[kk] == cl)
                     {
                        x = u.row.val[kk];
                        ll = kk;
                        break;
                     }
                  }

                  for(--kk; kk > m; --kk)
                  {
                     if(l_maxabs < spxAbs(u.row.val[kk]))
                        l_maxabs = spxAbs(u.row.val[kk]);
                  }

                  temp.s_max[k] = l_maxabs;
               }

               l_maxabs *= threshold;

               if(spxAbs(x) > l_maxabs)
               {
                  mkwtz = j;
                  rw = k;
                  ii = ll;

                  if(j <= count + 1)
                     break;
               }
            }
         }
      }
      else
      {
         ++count;
         continue;
      }

      assert(cl >= 0);

      removeDR(temp.pivot_col[cl]);
      initDR(temp.pivot_col[cl]);

      if(ii >= 0)
      {
         /*  Initialize selected pivot element
          */
         CLUFactor::Pring* pr;
         temp.pivot_row[rw].pos = ii - u.row.start[rw];
         temp.pivot_row[rw].mkwtz = mkwtz = (mkwtz - 1) * (count - 1);
         // ??? mkwtz originally was long,
         // maybe to avoid an overflow in this instruction?

         for(pr = temp.pivots.next; pr->idx >= 0; pr = pr->next)
         {
            if(pr->idx == rw || pr->mkwtz >= mkwtz)
               break;
         }

         pr = pr->prev;

         if(pr->idx != rw)
         {
            removeDR(temp.pivot_row[rw]);
            init2DR(temp.pivot_row[rw], *pr);
         }

         num++;

         if(num >= candidates)
            break;
      }
   }

   /*
    *     while(temp.temp.next->mkwtz < temp.temp.prev->mkwtz)
    *     {
    *     Pring   *pr;
    *     pr = temp.temp.prev;
    *     removeDR(*pr);
    *     init2DR (*pr, rowNZ[u.row.len[pr->idx]]);
    }
    */

   assert(row.perm[rw] < 0);

   assert(col.perm[cl] < 0);
}


/*
 *      Perform L and update loop for row r
 */
int CLUFactor::updateRow(int r,
                         int lv,
                         int prow,
                         int pcol,
                         Real pval,
                         Real eps)
{
   int fill;
   Real x, lx;
   int c, i, j, k, ll, m, n;

   n = u.row.start[r];
   m = --(u.row.len[r]) + n;

   /*  compute L vector entry and
    *  and remove pivot column form row file
    */

   for(j = m; u.row.idx[j] != pcol; --j)
      ;

   lx = u.row.val[j] / pval;

   l.val[lv] = lx;

   l.idx[lv] = r;

   ++lv;

   u.row.idx[j] = u.row.idx[m];

   u.row.val[j] = u.row.val[m];


   /*  update loop (I) and
    *  computing expected fill
    */
   fill = u.row.len[prow];

   for(j = m - 1; j >= n; --j)
   {
      c = u.row.idx[j];

      if(temp.s_mark[c])
      {
         /*  count fill elements.
          */
         temp.s_mark[c] = 0;
         --fill;

         /*  update row values
          */
         x = u.row.val[j] -= work[c] * lx;

         if(isZero(x, eps))
         {
            /* Eliminate zero from row r
             */
            --u.row.len[r];
            --m;
            u.row.val[j] = u.row.val[m];
            u.row.idx[j] = u.row.idx[m];

            /* Eliminate zero from column c
             */
            --(temp.s_cact[c]);
            k = --(u.col.len[c]) + u.col.start[c];

            for(i = k; u.col.idx[i] != r; --i)
               ;

            u.col.idx[i] = u.col.idx[k];
         }
      }
   }


   /*  create space for fill in row file
    */
   ll = u.row.len[r];

   if(ll + fill > u.row.max[r])
      remaxRow(r, ll + fill);

   ll += u.row.start[r];

   /*  fill creating update loop (II)
    */
   for(j = u.row.start[prow], m = j + u.row.len[prow]; j < m; ++j)
   {
      c = u.row.idx[j];

      if(temp.s_mark[c])
      {
         x = - work[c] * lx;

         if(isNotZero(x, eps))
         {
            /* produce fill element in row r
             */
            u.row.val[ll] = x;
            u.row.idx[ll] = c;
            ll++;
            u.row.len[r]++;

            /* produce fill element in column c
             */

            if(u.col.len[c] >= u.col.max[c])
               remaxCol(c, u.col.len[c] + 1);

            u.col.idx[u.col.start[c] + (u.col.len[c])++] = r;

            temp.s_cact[c]++;
         }
      }
      else
         temp.s_mark[c] = 1;
   }

   /*  move row to appropriate list.
    */
   removeDR(temp.pivot_row[r]);

   init2DR(temp.pivot_row[r], temp.pivot_rowNZ[u.row.len[r]]);

   assert(row.perm[r] < 0);

   temp.s_max[r] = -1;

   return lv;
}

/*
 *      Eliminate pivot element
 */
void CLUFactor::eliminatePivot(int prow, int pos, Real eps)
{
   int i, j, k, m = -1;
   int lv = -1;  // This value should never be used.
   int pcol;
   Real pval;
   int pbeg = u.row.start[prow];
   int plen = --(u.row.len[prow]);
   int pend = pbeg + plen;


   /*  extract pivot element   */
   i = pbeg + pos;
   pcol = u.row.idx[i];
   pval = u.row.val[i];
   removeDR(temp.pivot_col[pcol]);
   initDR(temp.pivot_col[pcol]);

   /*  remove pivot from pivot row     */
   u.row.idx[i] = u.row.idx[pend];
   u.row.val[i] = u.row.val[pend];

   /*  set pivot element and construct L vector */
   setPivot(temp.stage++, pcol, prow, pval);

   /**@todo If this test failes, lv has no value. I suppose that in this
    *       case none of the loops below that uses lv is executed.
    *       But this is unproven.
    */

   if(temp.s_cact[pcol] - 1 > 0)
      lv = makeLvec(temp.s_cact[pcol] - 1, prow);

   /*  init working vector,
    *  remove pivot row from working matrix
    *  and remove columns from list.
    */
   for(i = pbeg; i < pend; ++i)
   {
      j = u.row.idx[i];
      temp.s_mark[j] = 1;
      work[j] = u.row.val[i];
      removeDR(temp.pivot_col[j]);
      m = u.col.start[j] + u.col.len[j] - temp.s_cact[j];

      for(k = m; u.col.idx[k] != prow; ++k)
         ;

      u.col.idx[k] = u.col.idx[m];

      u.col.idx[m] = prow;

      temp.s_cact[j]--;
   }

   /*  perform L and update loop
    */
   for(i = u.col.len[pcol] - temp.s_cact[pcol];
         (m = u.col.idx[u.col.start[pcol] + i]) != prow;
         ++i)
   {
      assert(row.perm[m] < 0);
      assert(lv >= 0);
      updateRow(m, lv++, prow, pcol, pval, eps);
   }

   /*  skip pivot row  */

   m = u.col.len[pcol];

   for(++i; i < m; ++i)
   {
      assert(lv >= 0);
      updateRow(u.col.idx[u.col.start[pcol] + i], lv++, prow, pcol, pval, eps);
   }

   /*  remove pivot column from column file.
    */
   u.col.len[pcol] -= temp.s_cact[pcol];

   /*  clear working vector and reinsert columns to lists
    */
   for(i = u.row.start[prow], pend = i + plen; i < pend; ++i)
   {
      j = u.row.idx[i];
      work[j] = 0;
      temp.s_mark[j] = 0;
      init2DR(temp.pivot_col[j], temp.pivot_colNZ[temp.s_cact[j]]);
      assert(col.perm[j] < 0);
   }
}


/*
 *      Factorize nucleus.
 */
void CLUFactor::eliminateNucleus(const Real eps,
                                 const Real threshold)
{
   int r, c;
   CLUFactor::Pring* pivot;

   if(stat == SLinSolver::SINGULAR)
      return;

   temp.pivots.mkwtz = -1;

   temp.pivots.idx = -1;

   temp.pivots.pos = -1;

   while(temp.stage < thedim - 1)
   {
#ifndef NDEBUG
      int i;
      // CLUFactorIsConsistent(fac);

      for(i = 0; i < thedim; ++i)
         if(col.perm[i] < 0)
            assert(temp.s_mark[i] == 0);

#endif

      if(temp.pivot_rowNZ[1].next != &(temp.pivot_rowNZ[1]))
         /* row singleton available */
         eliminateRowSingletons();
      else if(temp.pivot_colNZ[1].next != &(temp.pivot_colNZ[1]))
         /* column singleton available */
         eliminateColSingletons();
      else
      {
         initDR(temp.pivots);
         selectPivots(threshold);

         assert(temp.pivots.next != &temp.pivots &&
                "ERROR: no pivot element selected");

         for(pivot = temp.pivots.next; pivot != &temp.pivots;
               pivot = pivot->next)
         {
            eliminatePivot(pivot->idx, pivot->pos, eps);
         }
      }

      if(temp.pivot_rowNZ->next != temp.pivot_rowNZ ||
            temp.pivot_colNZ->next != temp.pivot_colNZ)
      {
         stat = SLinSolver::SINGULAR;
         return;
      }
   }

   if(temp.stage < thedim)
   {
      /*      Eliminate remaining element.
       *      Note, that this must be both, column and row singleton.
       */
      assert(temp.pivot_rowNZ[1].next != &(temp.pivot_rowNZ[1]) &&
             "ERROR: one row must be left");
      assert(temp.pivot_colNZ[1].next != &(temp.pivot_colNZ[1]) &&
             "ERROR: one col must be left");
      r = temp.pivot_rowNZ[1].next->idx;
      c = temp.pivot_colNZ[1].next->idx;
      u.row.len[r] = 0;
      u.col.len[c]--;
      setPivot(temp.stage, c, r, u.row.val[u.row.start[r]]);
   }
}

/*****************************************************************************/

int CLUFactor::setupColVals()
{
   int i;
   int n = thedim;

   if(u.col.val != 0)
      spx_free(u.col.val);

   // if we would know the old size of u.col.val, this could be a realloc.
   spx_alloc(u.col.val, u.col.size);

   for(i = 0; i < thedim; i++)
      u.col.len[i] = 0;

   maxabs = REAL(0.0);

   for(i = 0; i < thedim; i++)
   {
      int     k   = u.row.start[i];
      int*    idx = &u.row.idx[k];
      Real*   val = &u.row.val[k];
      int     len = u.row.len[i];

      n += len;

      while(len-- > 0)
      {
         assert((*idx >= 0) && (*idx < thedim));

         k = u.col.start[*idx] + u.col.len[*idx];

         assert((k >= 0) && (k < u.col.size));

         u.col.len[*idx]++;

         assert(u.col.len[*idx] <= u.col.max[*idx]);

         u.col.idx[k] = i;
         u.col.val[k] = *val;

         if(spxAbs(*val) > maxabs)
            maxabs = spxAbs(*val);

         idx++;

         val++;
      }
   }

   return n;
}

/*****************************************************************************/

#ifdef WITH_L_ROWS
void CLUFactor::setupRowVals()
{
   int   i, j, k, m;
   int   vecs, mem;
   int*  l_row;
   int*  idx;
   Real* val;
   int*  beg;
   int*  l_ridx;
   Real* l_rval;
   int*  l_rbeg;
   int*  rorig;
   int*  rrorig;
   int*  rperm;
   int*  rrperm;

   vecs  = l.firstUpdate;
   l_row = l.row;
   idx   = l.idx;
   val   = l.val;
   beg   = l.start;
   mem   = beg[vecs];

   if(l.rval)
      spx_free(l.rval);

   if(l.ridx)
      spx_free(l.ridx);

   if(l.rbeg)
      spx_free(l.rbeg);

   if(l.rorig)
      spx_free(l.rorig);

   if(l.rperm)
      spx_free(l.rperm);

   spx_alloc(l.rval, mem);

   spx_alloc(l.ridx, mem);

   spx_alloc(l.rbeg, thedim + 1);

   spx_alloc(l.rorig, thedim);

   spx_alloc(l.rperm, thedim);

   l_ridx = l.ridx;

   l_rval = l.rval;

   l_rbeg = l.rbeg;

   rorig  = l.rorig;

   rrorig = row.orig;

   rperm  = l.rperm;

   rrperm = row.perm;

   for(i = thedim; i--; *l_rbeg++ = 0)
   {
      *rorig++ = *rrorig++;
      *rperm++ = *rrperm++;
   }

   *l_rbeg = 0;

   l_rbeg = l.rbeg + 1;

   for(i = mem; i--;)
      l_rbeg[*idx++]++;

   idx = l.idx;

   for(m = 0, i = thedim; i--; l_rbeg++)
   {
      j = *l_rbeg;
      *l_rbeg = m;
      m += j;
   }

   assert(m == mem);

   l_rbeg = l.rbeg + 1;

   for(i = j = 0; i < vecs; ++i)
   {
      m = l_row[i];
      assert(idx == &l.idx[l.start[i]]);

      for(; j < beg[i + 1]; j++)
      {
         k = l_rbeg[*idx++]++;
         assert(k < mem);
         l_ridx[k] = m;
         l_rval[k] = *val++;
      }
   }

   assert(l.rbeg[thedim] == mem);

   assert(l.rbeg[0] == 0);
}

#endif

/*****************************************************************************/

void CLUFactor::factor(const SVector** vec,          ///< Array of column vector pointers
                       Real            threshold,    ///< pivoting threshold
                       Real            eps)          ///< epsilon for zero detection
{

   factorTime->start();

   stat = SLinSolver::OK;

   l.start[0]    = 0;
   l.firstUpdate = 0;
   l.firstUnused = 0;

   temp.init(thedim);
   initPerm();

   initFactorMatrix(vec, eps);

   if(stat)
      goto TERMINATE;

   //   initMaxabs = initMaxabs;

   colSingletons();

   if(stat != SLinSolver::OK)
      goto TERMINATE;

   rowSingletons();

   if(stat != SLinSolver::OK)
      goto TERMINATE;

   if(temp.stage < thedim)
   {
      initFactorRings();
      eliminateNucleus(eps, threshold);
      freeFactorRings();
   }

TERMINATE:

   l.firstUpdate = l.firstUnused;

   if(stat == SLinSolver::OK)
   {
#ifdef WITH_L_ROWS
      setupRowVals();
#endif
      nzCnt = setupColVals();
   }

   factorTime->stop();

   factorCount++;
}

void CLUFactor::dump() const
{
   int i, j, k;

   // Dump regardless of the verbosity level if this method is called;

   /*  Dump U:
    */

   for(i = 0; i < thedim; ++i)
   {
      if(row.perm[i] >= 0)
         std::cout << "DCLUFA01 diag[" << i << "]: [" << col.orig[row.perm[i]]
                   << "] = " << diag[i] << std::endl;

      for(j = 0; j < u.row.len[i]; ++j)
         std::cout << "DCLUFA02   u[" << i << "]: ["
                   << u.row.idx[u.row.start[i] + j] << "] = "
                   << u.row.val[u.row.start[i] + j] << std::endl;
   }

   /*  Dump L:
    */
   for(i = 0; i < thedim; ++i)
   {
      for(j = 0; j < l.firstUnused; ++j)
         if(col.orig[row.perm[l.row[j]]] == i)
         {
            std::cout << "DCLUFA03 l[" << i << "]" << std::endl;

            for(k = l.start[j]; k < l.start[j + 1]; ++k)
               std::cout << "DCLUFA04   l[" << k - l.start[j]
                         << "]:  [" << l.idx[k]
                         << "] = "  << l.val[k] << std::endl;

            break;
         }
   }

   return;
}

/*****************************************************************************/
/*
 *      Ensure that row memory is at least size.
 */
void CLUFactor::minRowMem(int size)
{

   if(u.row.size < size)
   {
      u.row.size = size;
      spx_realloc(u.row.val, size);
      spx_realloc(u.row.idx, size);
   }
}

/*****************************************************************************/
/*
 *      Ensure that column memory is at least size.
 */
void CLUFactor::minColMem(int size)
{

   if(u.col.size < size)
   {
      u.col.size = size;
      spx_realloc(u.col.idx, size);
   }
}

void CLUFactor::forestMinColMem(int size)
{

   if(u.col.size < size)
   {
      u.col.size = size;
      spx_realloc(u.col.idx, size);
      spx_realloc(u.col.val, size);
   }
}

void CLUFactor::minLMem(int size)
{

   if(size > l.size)
   {
      l.size = int(0.2 * l.size + size);
      spx_realloc(l.val, l.size);
      spx_realloc(l.idx, l.size);
   }
}


int CLUFactor::makeLvec(int p_len, int p_row)
{

   if(l.firstUnused >= l.startSize)
   {
      l.startSize += 100;
      spx_realloc(l.start, l.startSize);
   }

   int* p_lrow = l.row;

   int* p_lbeg = l.start;
   int first   = p_lbeg[l.firstUnused];

   assert(p_len > 0 && "ERROR: no empty columns allowed in L vectors");

   minLMem(first + p_len);
   p_lrow[l.firstUnused] = p_row;
   l.start[++(l.firstUnused)] = first + p_len;

   assert(l.start[l.firstUnused] <= l.size);
   assert(l.firstUnused <= l.startSize);
   return first;
}


/*****************************************************************************/

bool CLUFactor::isConsistent() const
{
#ifdef ENABLE_CONSISTENCY_CHECKS
   int              i, j, k, ll;
   Dring*            ring;
   CLUFactor::Pring* pring;

   /*  Consistency only relevant for real factorizations
    */

   if(stat)
      return true;

   /*  Test column ring list consistency.
    */
   i = 0;

   for(ring = u.col.list.next; ring != &(u.col.list); ring = ring->next)
   {
      assert(ring->idx >= 0);
      assert(ring->idx < thedim);
      assert(ring->prev->next == ring);

      if(ring != u.col.list.next)
      {
         assert(u.col.start[ring->prev->idx] + u.col.max[ring->prev->idx]
                == u.col.start[ring->idx]);
      }

      ++i;
   }

   assert(i == thedim);

   assert(u.col.start[ring->prev->idx] + u.col.max[ring->prev->idx]
          == u.col.used);


   /*  Test row ring list consistency.
    */
   i = 0;

   for(ring = u.row.list.next; ring != &(u.row.list); ring = ring->next)
   {
      assert(ring->idx >= 0);
      assert(ring->idx < thedim);
      assert(ring->prev->next == ring);
      assert(u.row.start[ring->prev->idx] + u.row.max[ring->prev->idx]
             == u.row.start[ring->idx]);
      ++i;
   }

   assert(i == thedim);

   assert(u.row.start[ring->prev->idx] + u.row.max[ring->prev->idx]
          == u.row.used);


   /*  Test consistency of individual svectors.
    */

   for(i = 0; i < thedim; ++i)
   {
      assert(u.row.max[i] >= u.row.len[i]);
      assert(u.col.max[i] >= u.col.len[i]);
   }


   /*  Test consistency of column file to row file of U
    */
   for(i = 0; i < thedim; ++i)
   {
      for(j = u.row.start[i] + u.row.len[i] - 1; j >= u.row.start[i]; j--)
      {
         k = u.row.idx[j];

         for(ll = u.col.start[k] + u.col.len[k] - 1; ll >= u.col.start[k]; ll--)
         {
            if(u.col.idx[ll] == i)
               break;
         }

         assert(u.col.idx[ll] == i);

         if(row.perm[i] < 0)
         {
            assert(col.perm[k] < 0);
         }
         else
         {
            assert(col.perm[k] < 0 || col.perm[k] > row.perm[i]);
         }
      }
   }

   /*  Test consistency of row file to column file of U
    */
   for(i = 0; i < thedim; ++i)
   {
      for(j = u.col.start[i] + u.col.len[i] - 1; j >= u.col.start[i]; j--)
      {
         k = u.col.idx[j];

         for(ll = u.row.start[k] + u.row.len[k] - 1; ll >= u.row.start[k]; ll--)
         {
            if(u.row.idx[ll] == i)
               break;
         }

         assert(u.row.idx[ll] == i);

         assert(col.perm[i] < 0 || row.perm[k] < col.perm[i]);
      }
   }

   /*  Test consistency of nonzero count lists
    */
   if(temp.pivot_colNZ && temp.pivot_rowNZ)
   {
      for(i = 0; i < thedim - temp.stage; ++i)
      {
         for(pring = temp.pivot_rowNZ[i].next; pring != &(temp.pivot_rowNZ[i]); pring = pring->next)
         {
            assert(row.perm[pring->idx] < 0);
         }

         for(pring = temp.pivot_colNZ[i].next; pring != &(temp.pivot_colNZ[i]); pring = pring->next)
         {
            assert(col.perm[pring->idx] < 0);
         }
      }
   }

#endif // CONSISTENCY_CHECKS

   return true;
}

void CLUFactor::solveUright(Real* wrk, Real* vec) const
{

   for(int i = thedim - 1; i >= 0; i--)
   {
      int  r = row.orig[i];
      int  c = col.orig[i];
      Real x = wrk[c] = diag[r] * vec[r];

      vec[r] = 0.0;

      if(x != 0.0)
         //if (isNotZero(x))
      {
         for(int j = u.col.start[c]; j < u.col.start[c] + u.col.len[c]; j++)
            vec[u.col.idx[j]] -= x * u.col.val[j];
      }
   }
}

int CLUFactor::solveUrightEps(Real* vec, int* nonz, Real eps, Real* rhs)
{
   int i, j, r, c, n;
   int* rorig, *corig;
   int* cidx, *clen, *cbeg;
   Real* cval;
   Real x;

   int* idx;
   Real* val;

   rorig = row.orig;
   corig = col.orig;

   cidx = u.col.idx;
   cval = u.col.val;
   clen = u.col.len;
   cbeg = u.col.start;

   n = 0;

   for(i = thedim - 1; i >= 0; --i)
   {
      r = rorig[i];
      x = diag[r] * rhs[r];

      if(isNotZero(x, eps))
      {
         c = corig[i];
         vec[c] = x;
         nonz[n++] = c;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
            rhs[*idx++] -= x * (*val++);
      }
   }

   return n;
}

void CLUFactor::solveUright2(
   Real* p_work1, Real* vec1, Real* p_work2, Real* vec2)
{
   int i, j, r, c;
   int* rorig, *corig;
   int* cidx, *clen, *cbeg;
   Real* cval;
   Real x1, x2;

   int* idx;
   Real* val;

   rorig = row.orig;
   corig = col.orig;

   cidx = u.col.idx;
   cval = u.col.val;
   clen = u.col.len;
   cbeg = u.col.start;

   for(i = thedim - 1; i >= 0; --i)
   {
      r = rorig[i];
      c = corig[i];
      p_work1[c] = x1 = diag[r] * vec1[r];
      p_work2[c] = x2 = diag[r] * vec2[r];
      vec1[r] = vec2[r] = 0;

      if(x1 != 0.0 && x2 != 0.0)
      {
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
         {
            vec1[*idx] -= x1 * (*val);
            vec2[*idx++] -= x2 * (*val++);
         }
      }
      else if(x1 != 0.0)
      {
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
            vec1[*idx++] -= x1 * (*val++);
      }
      else if(x2 != 0.0)
      {
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
            vec2[*idx++] -= x2 * (*val++);
      }
   }
}

int CLUFactor::solveUright2eps(
   Real* p_work1, Real* vec1, Real* p_work2, Real* vec2,
   int* nonz, Real eps)
{
   int i, j, r, c, n;
   int* rorig, *corig;
   int* cidx, *clen, *cbeg;
   bool notzero1, notzero2;
   Real* cval;
   Real x1, x2;

   int* idx;
   Real* val;

   rorig = row.orig;
   corig = col.orig;

   cidx = u.col.idx;
   cval = u.col.val;
   clen = u.col.len;
   cbeg = u.col.start;

   n = 0;

   for(i = thedim - 1; i >= 0; --i)
   {
      c = corig[i];
      r = rorig[i];
      p_work1[c] = x1 = diag[r] * vec1[r];
      p_work2[c] = x2 = diag[r] * vec2[r];
      vec1[r] = vec2[r] = 0;
      notzero1 = isNotZero(x1, eps);
      notzero2 = isNotZero(x2, eps);

      if(notzero1 && notzero2)
      {
         *nonz++ = c;
         n++;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
         {
            vec1[*idx] -= x1 * (*val);
            vec2[*idx++] -= x2 * (*val++);
         }
      }
      else if(notzero1)
      {
         p_work2[c] = 0.0;
         *nonz++ = c;
         n++;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
            vec1[*idx++] -= x1 * (*val++);
      }
      else if(notzero2)
      {
         p_work1[c] = 0.0;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
            vec2[*idx++] -= x2 * (*val++);
      }
      else
      {
         p_work1[c] = 0.0;
         p_work2[c] = 0.0;
      }
   }

   return n;
}

void CLUFactor::solveLright(Real* vec)
{
   int i, j, k;
   int end;
   Real x;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   for(i = 0; i < end; ++i)
   {
      if((x = vec[lrow[i]]) != 0.0)
      {
         MSG_DEBUG(std::cout << "y" << lrow[i] << "=" << vec[lrow[i]] << std::endl;)

         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            MSG_DEBUG(std::cout << "                         -> y" << *idx << " -= " << x << " * " << *val <<
                      " = " << x * (*val) << "    -> " << vec[*idx] - x * (*val) << std::endl;)
            vec[*idx++] -= x * (*val++);
         }
      }
   }

   if(l.updateType)                      /* Forest-Tomlin Updates */
   {
      MSG_DEBUG(std::cout << "performing FT updates..." << std::endl;)

      end = l.firstUnused;

      for(; i < end; ++i)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         StableSum<Real> tmp(-vec[lrow[i]]);

         for(j = lbeg[i + 1]; j > k; --j)
            tmp += vec[*idx++] * (*val++);

         vec[lrow[i]] = -tmp;

         MSG_DEBUG(std::cout << "y" << lrow[i] << "=" << vec[lrow[i]] << std::endl;)
      }

      MSG_DEBUG(std::cout << "finished FT updates." << std::endl;)
   }
}

void CLUFactor::solveLright2(Real* vec1, Real* vec2)
{
   int i, j, k;
   int end;
   Real x2;
   Real x1;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   for(i = 0; i < end; ++i)
   {
      x1 = vec1[lrow[i]];
      x2 = vec2[lrow[i]];

      if(x1 != 0 && x2 != 0.0)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            vec1[*idx] -= x1 * (*val);
            vec2[*idx++] -= x2 * (*val++);
         }
      }
      else if(x1 != 0.0)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
            vec1[*idx++] -= x1 * (*val++);
      }
      else if(x2 != 0.0)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
            vec2[*idx++] -= x2 * (*val++);
      }
   }

   if(l.updateType)                      /* Forest-Tomlin Updates */
   {
      end = l.firstUnused;

      for(; i < end; ++i)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         StableSum<Real> tmp1(-vec1[lrow[i]]);
         StableSum<Real> tmp2(-vec2[lrow[i]]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            tmp1 += vec1[*idx] * (*val);
            tmp2 += vec2[*idx++] * (*val++);
         }

         vec1[lrow[i]] = -tmp1;
         vec2[lrow[i]] = -tmp2;
      }
   }
}

void CLUFactor::solveUpdateRight(Real* vec)
{
   int i, j, k;
   int end;
   Real x;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   assert(!l.updateType);               /* no Forest-Tomlin Updates */

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUnused;

   for(i = l.firstUpdate; i < end; ++i)
   {
      if((x = vec[lrow[i]]) != 0.0)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
            vec[*idx++] -= x * (*val++);
      }
   }
}

void CLUFactor::solveUpdateRight2(Real* vec1, Real* vec2)
{
   int i, j, k;
   int end;
   Real x1, x2;
   Real* lval;
   int* lrow, *lidx;
   int* lbeg;

   int* idx;
   Real* val;

   assert(!l.updateType);               /* no Forest-Tomlin Updates */

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUnused;

   for(i = l.firstUpdate; i < end; ++i)
   {
      x1 = vec1[lrow[i]];
      x2 = vec2[lrow[i]];

      if(x1 != 0.0 && x2 != 0.0)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            vec1[*idx] -= x1 * (*val);
            vec2[*idx++] -= x2 * (*val++);
         }
      }
      else if(x1 != 0.0)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
            vec1[*idx++] -= x1 * (*val++);
      }
      else if(x2 != 0.0)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
            vec2[*idx++] -= x2 * (*val++);
      }
   }
}

int CLUFactor::solveRight4update(Real* vec, int* nonz, Real eps,
                                 Real* rhs, Real* forest, int* forestNum, int* forestIdx)
{
   solveLright(rhs);

   if(forest)
   {
      int n = 0;

      for(int i = 0; i < thedim; i++)
      {
         forestIdx[n] = i;
         forest[i]    = rhs[i];
         n           += rhs[i] != 0.0 ? 1 : 0;
      }

      *forestNum = n;
   }

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      solveUright(vec, rhs);
      solveUpdateRight(vec);
      return 0;
   }
   else
      return solveUrightEps(vec, nonz, eps, rhs);
}

void CLUFactor::solveRight(Real* vec, Real* rhs)
{
   solveLright(rhs);
   solveUright(vec, rhs);

   if(!l.updateType)             /* no Forest-Tomlin Updates */
      solveUpdateRight(vec);
}

int CLUFactor::solveRight2update(Real* vec1,
                                 Real* vec2,
                                 Real* rhs1,
                                 Real* rhs2,
                                 int* nonz,
                                 Real eps,
                                 Real* forest,
                                 int* forestNum,
                                 int* forestIdx)
{
   solveLright2(rhs1, rhs2);

   if(forest)
   {
      int n = 0;

      for(int i = 0; i < thedim; i++)
      {
         forestIdx[n] = i;
         forest[i]    = rhs1[i];
         n           += rhs1[i] != 0.0 ? 1 : 0;
      }

      *forestNum = n;
   }

   if(!l.updateType)            /* no Forest-Tomlin Updates */
   {
      solveUright2(vec1, rhs1, vec2, rhs2);
      solveUpdateRight2(vec1, vec2);
      return 0;
   }
   else
      return solveUright2eps(vec1, rhs1, vec2, rhs2, nonz, eps);
}

void CLUFactor::solveRight2(
   Real* vec1,
   Real* vec2,
   Real* rhs1,
   Real* rhs2)
{
   solveLright2(rhs1, rhs2);

   if(l.updateType)              /* Forest-Tomlin Updates */
      solveUright2(vec1, rhs1, vec2, rhs2);
   else
   {
      solveUright2(vec1, rhs1, vec2, rhs2);
      solveUpdateRight2(vec1, vec2);
   }
}

/*****************************************************************************/
void CLUFactor::solveUleft(Real* p_work, Real* vec)
{
   for(int i = 0; i < thedim; ++i)
   {
      int  c  = col.orig[i];
      int  r  = row.orig[i];

      assert(c >= 0);    // Inna/Tobi: otherwise, vec[c] would be strange...
      assert(r >= 0);    // Inna/Tobi: otherwise, diag[r] would be strange...

      Real x  = vec[c];

      ASSERT_WARN("WSOLVE01", spxAbs(vec[c]) < 1e40);

      vec[c]  = 0.0;

      if(x != 0.0)
      {
         ASSERT_WARN("WSOLVE02", spxAbs(diag[r]) < 1e40);

         x        *= diag[r];
         p_work[r] = x;

         int end = u.row.start[r] + u.row.len[r];

         for(int m = u.row.start[r]; m < end; m++)
         {
            vec[u.row.idx[m]] -= x * u.row.val[m];
            ASSERT_WARN("WSOLVE03", spxAbs(vec[u.row.idx[m]]) < 1e40);
         }
      }
   }
}

void CLUFactor::solveUleft2(
   Real* p_work1, Real* vec1, Real* p_work2, Real* vec2)
{
   Real x1;
   Real x2;
   int i, k, r, c;
   int* rorig, *corig;
   int* ridx, *rlen, *rbeg, *idx;
   Real* rval, *val;

   rorig = row.orig;
   corig = col.orig;

   ridx = u.row.idx;
   rval = u.row.val;
   rlen = u.row.len;
   rbeg = u.row.start;

   for(i = 0; i < thedim; ++i)
   {
      c = corig[i];
      r = rorig[i];
      x1 = vec1[c];
      x2 = vec2[c];

      if((x1 != 0.0) && (x2 != 0.0))
      {
         x1 *= diag[r];
         x2 *= diag[r];
         p_work1[r] = x1;
         p_work2[r] = x2;
         k = rbeg[r];
         idx = &ridx[k];
         val = &rval[k];

         for(int m = rlen[r]; m != 0; --m)
         {
            vec1[*idx] -= x1 * (*val);
            vec2[*idx] -= x2 * (*val++);
            idx++;
         }
      }
      else if(x1 != 0.0)
      {
         x1 *= diag[r];
         p_work1[r] = x1;
         k = rbeg[r];
         idx = &ridx[k];
         val = &rval[k];

         for(int m = rlen[r]; m != 0; --m)
            vec1[*idx++] -= x1 * (*val++);
      }
      else if(x2 != 0.0)
      {
         x2 *= diag[r];
         p_work2[r] = x2;
         k = rbeg[r];
         idx = &ridx[k];
         val = &rval[k];

         for(int m = rlen[r]; m != 0; --m)
            vec2[*idx++] -= x2 * (*val++);
      }
   }
}

int CLUFactor::solveLleft2forest(
   Real* vec1,
   int* /* nonz */,
   Real* vec2,
   Real /* eps */)
{
   int i;
   int j;
   int k;
   int end;
   Real x1, x2;
   Real* lval, *val;
   int* lidx, *idx, *lrow;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   for(i = l.firstUnused - 1; i >= end; --i)
   {
      j = lrow[i];
      x1 = vec1[j];
      x2 = vec2[j];

      if(x1 != 0.0)
      {
         if(x2 != 0.0)
         {
            k = lbeg[i];
            val = &lval[k];
            idx = &lidx[k];

            for(j = lbeg[i + 1]; j > k; --j)
            {
               vec1[*idx] -= x1 * (*val);
               vec2[*idx++] -= x2 * (*val++);
            }
         }
         else
         {
            k = lbeg[i];
            val = &lval[k];
            idx = &lidx[k];

            for(j = lbeg[i + 1]; j > k; --j)
               vec1[*idx++] -= x1 * (*val++);
         }
      }
      else if(x2 != 0.0)
      {
         k = lbeg[i];
         val = &lval[k];
         idx = &lidx[k];

         for(j = lbeg[i + 1]; j > k; --j)
            vec2[*idx++] -= x2 * (*val++);
      }
   }

   return 0;
}

void CLUFactor::solveLleft2(
   Real* vec1,
   int* /* nonz */,
   Real* vec2,
   Real /* eps */)
{
   int i, j, k, r;
   int x1not0, x2not0;
   Real x1, x2;

   Real* rval, *val;
   int* ridx, *idx;
   int* rbeg;
   int* rorig;

   ridx  = l.ridx;
   rval  = l.rval;
   rbeg  = l.rbeg;
   rorig = l.rorig;

#ifndef WITH_L_ROWS
   Real*   lval  = l.val;
   int*    lidx  = l.idx;
   int*    lrow  = l.row;
   int*    lbeg  = l.start;

   i = l.firstUpdate - 1;

   for(; i >= 0; --i)
   {
      k = lbeg[i];
      val = &lval[k];
      idx = &lidx[k];
      x1 = 0;
      x2 = 0;

      for(j = lbeg[i + 1]; j > k; --j)
      {
         x1 += vec1[*idx] * (*val);
         x2 += vec2[*idx++] * (*val++);
      }

      vec1[lrow[i]] -= x1;

      vec2[lrow[i]] -= x2;
   }

#else

   for(i = thedim; i--;)
   {
      r = rorig[i];
      x1 = vec1[r];
      x2 = vec2[r];
      x1not0 = (x1 != 0);
      x2not0 = (x2 != 0);

      if(x1not0 && x2not0)
      {
         k = rbeg[r];
         j = rbeg[r + 1] - k;
         val = &rval[k];
         idx = &ridx[k];

         while(j-- > 0)
         {
            assert(row.perm[*idx] < i);
            vec1[*idx] -= x1 * *val;
            vec2[*idx++] -= x2 * *val++;
         }
      }
      else if(x1not0)
      {
         k = rbeg[r];
         j = rbeg[r + 1] - k;
         val = &rval[k];
         idx = &ridx[k];

         while(j-- > 0)
         {
            assert(row.perm[*idx] < i);
            vec1[*idx++] -= x1 * *val++;
         }
      }
      else if(x2not0)
      {
         k = rbeg[r];
         j = rbeg[r + 1] - k;
         val = &rval[k];
         idx = &ridx[k];

         while(j-- > 0)
         {
            assert(row.perm[*idx] < i);
            vec2[*idx++] -= x2 * *val++;
         }
      }
   }

#endif
}

int CLUFactor::solveLleftForest(Real* vec, int* /* nonz */, Real /* eps */)
{
   int i, j, k, end;
   Real x;
   Real* val, *lval;
   int* idx, *lidx, *lrow, *lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   for(i = l.firstUnused - 1; i >= end; --i)
   {
      if((x = vec[lrow[i]]) != 0.0)
      {
         k = lbeg[i];
         val = &lval[k];
         idx = &lidx[k];

         for(j = lbeg[i + 1]; j > k; --j)
            vec[*idx++] -= x * (*val++);
      }
   }

   return 0;
}

void CLUFactor::solveLleft(Real* vec) const
{

#ifndef WITH_L_ROWS
   int*  idx;
   Real* val;
   Real* lval  = l.val;
   int*  lidx  = l.idx;
   int*  lrow  = l.row;
   int*  lbeg  = l.start;

   for(int i = l.firstUpdate - 1; i >= 0; --i)
   {
      int k = lbeg[i];
      val = &lval[k];
      idx = &lidx[k];
      x = 0;

      for(int j = lbeg[i + 1]; j > k; --j)
         x += vec[*idx++] * (*val++);

      vec[lrow[i]] -= x;
   }

#else

   for(int i = thedim - 1; i >= 0; --i)
   {
      int  r = l.rorig[i];
      Real x = vec[r];

      if(x != 0.0)
      {
         for(int k = l.rbeg[r]; k < l.rbeg[r + 1]; k++)
         {
            int j = l.ridx[k];

            assert(l.rperm[j] < i);

            vec[j] -= x * l.rval[k];
         }
      }
   }

#endif // WITH_L_ROWS
}

int CLUFactor::solveLleftEps(Real* vec, int* nonz, Real eps)
{
   int i, j, k, n;
   int r;
   Real x;
   Real* rval, *val;
   int* ridx, *idx;
   int* rbeg;
   int* rorig;

   ridx = l.ridx;
   rval = l.rval;
   rbeg = l.rbeg;
   rorig = l.rorig;
   n = 0;
#ifndef WITH_L_ROWS
   Real* lval = l.val;
   int*  lidx = l.idx;
   int*  lrow = l.row;
   int*  lbeg = l.start;

   for(i = l.firstUpdate - 1; i >= 0; --i)
   {
      k = lbeg[i];
      val = &lval[k];
      idx = &lidx[k];
      x = 0;

      for(j = lbeg[i + 1]; j > k; --j)
         x += vec[*idx++] * (*val++);

      vec[lrow[i]] -= x;
   }

#else

   for(i = thedim; i--;)
   {
      r = rorig[i];
      x = vec[r];

      if(isNotZero(x, eps))
      {
         *nonz++ = r;
         n++;
         k = rbeg[r];
         j = rbeg[r + 1] - k;
         val = &rval[k];
         idx = &ridx[k];

         while(j-- > 0)
         {
            assert(row.perm[*idx] < i);
            vec[*idx++] -= x * *val++;
         }
      }
      else
         vec[r] = 0;
   }

#endif

   return n;
}

void CLUFactor::solveUpdateLeft(Real* vec)
{
   int i, j, k, end;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   assert(!l.updateType);               /* Forest-Tomlin Updates */

   end = l.firstUpdate;

   for(i = l.firstUnused - 1; i >= end; --i)
   {
      k = lbeg[i];
      val = &lval[k];
      idx = &lidx[k];
      StableSum<Real> tmp(-vec[lrow[i]]);

      for(j = lbeg[i + 1]; j > k; --j)
         tmp += vec[*idx++] * (*val++);

      vec[lrow[i]] = -tmp;
   }
}

void CLUFactor::solveUpdateLeft2(Real* vec1, Real* vec2)
{
   int i, j, k, end;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   assert(!l.updateType);               /* Forest-Tomlin Updates */

   end = l.firstUpdate;

   for(i = l.firstUnused - 1; i >= end; --i)
   {
      k = lbeg[i];
      val = &lval[k];
      idx = &lidx[k];

      StableSum<Real> tmp1(-vec1[lrow[i]]);
      StableSum<Real> tmp2(-vec2[lrow[i]]);

      for(j = lbeg[i + 1]; j > k; --j)
      {
         tmp1 += vec1[*idx] * (*val);
         tmp2 += vec2[*idx++] * (*val++);
      }

      vec1[lrow[i]] = -tmp1;
      vec2[lrow[i]] = -tmp2;
   }
}

int CLUFactor::solveUpdateLeft(Real eps, Real* vec, int* nonz, int n)
{
   int i, j, k, end;
   Real y;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   assert(!l.updateType);               /* no Forest-Tomlin Updates! */

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   for(i = l.firstUnused - 1; i >= end; --i)
   {
      k = lbeg[i];
      assert(k >= 0 && k < l.size);
      val = &lval[k];
      idx = &lidx[k];

      k = lrow[i];

      y = vec[k];
      StableSum<Real> tmp(-y);

      for(j = lbeg[i + 1]; j > k; --j)
      {
         assert(*idx >= 0 && *idx < thedim);
         tmp += vec[*idx++] * (*val++);
      }

      if(y == 0)
      {
         y = -tmp;

         if(isNotZero(y, eps))
         {
            nonz[n++] = k;
            vec[k] = y;
         }
      }
      else
      {
         y = -tmp;
         vec[k] = (y != 0) ? y : SOPLEX_FACTOR_MARKER;
      }
   }

   return n;
}

void CLUFactor::solveLeft(Real* vec, Real* rhs)
{

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      solveUpdateLeft(rhs);
      solveUleft(vec, rhs);
      solveLleft(vec);
   }
   else
   {
      solveUleft(vec, rhs);
      //@todo is 0.0 the right value for eps here ?
      solveLleftForest(vec, 0, 0.0);
      solveLleft(vec);
   }
}

int CLUFactor::solveLeftEps(Real* vec, Real* rhs, int* nonz, Real eps)
{

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      solveUpdateLeft(rhs);
      solveUleft(vec, rhs);
      return solveLleftEps(vec, nonz, eps);
   }
   else
   {
      solveUleft(vec, rhs);
      solveLleftForest(vec, nonz, eps);
      return solveLleftEps(vec, nonz, eps);
   }
}

int CLUFactor::solveLeft2(
   Real* vec1,
   int* nonz,
   Real* vec2,
   Real eps,
   Real* rhs1,
   Real* rhs2)
{

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      solveUpdateLeft2(rhs1, rhs2);
      solveUleft2(vec1, rhs1, vec2, rhs2);
      solveLleft2(vec1, nonz, vec2, eps);
      return 0;
   }
   else
   {
      solveUleft2(vec1, rhs1, vec2, rhs2);
      solveLleft2forest(vec1, nonz, vec2, eps);
      solveLleft2(vec1, nonz, vec2, eps);
      return 0;
   }
}

int CLUFactor::solveUleft(Real eps,
                          Real* vec, int* vecidx,
                          Real* rhs, int* rhsidx, int rhsn)
{
   Real x, y;
   int i, j, k, n, r, c;
   int* rorig, *corig, *cperm;
   int* ridx, *rlen, *rbeg, *idx;
   Real* rval, *val;

   rorig = row.orig;
   corig = col.orig;
   cperm = col.perm;

   /*  move rhsidx to a heap
    */

   for(i = 0; i < rhsn;)
      enQueueMin(rhsidx, &i, cperm[rhsidx[i]]);

   ridx = u.row.idx;

   rval = u.row.val;

   rlen = u.row.len;

   rbeg = u.row.start;

   n = 0;

   while(rhsn > 0)
   {
      i = deQueueMin(rhsidx, &rhsn);
      assert(i >= 0 && i < thedim);
      c = corig[i];
      assert(c >= 0 && c < thedim);
      x = rhs[c];
      rhs[c] = 0;

      if(isNotZero(x, eps))
      {
         r = rorig[i];
         assert(r >= 0 && r < thedim);
         vecidx[n++] = r;
         x *= diag[r];
         vec[r] = x;
         k = rbeg[r];
         assert(k >= 0 && k < u.row.size);
         idx = &ridx[k];
         val = &rval[k];

         for(int m = rlen[r]; m; --m)
         {
            j = *idx++;
            assert(j >= 0 && j < thedim);
            y = rhs[j];

            if(y == 0)
            {
               y = -x * (*val++);

               if(isNotZero(y, eps))
               {
                  rhs[j] = y;
                  enQueueMin(rhsidx, &rhsn, cperm[j]);
               }
            }
            else
            {
               y -= x * (*val++);
               rhs[j] = (y != 0) ? y : SOPLEX_FACTOR_MARKER;
            }
         }
      }
   }

   return n;
}


void CLUFactor::solveUleftNoNZ(Real eps, Real* vec,
                               Real* rhs, int* rhsidx, int rhsn)
{
   Real x, y;
   int i, j, k, r, c;
   int* rorig, *corig, *cperm;
   int* ridx, *rlen, *rbeg, *idx;
   Real* rval, *val;

   rorig = row.orig;
   corig = col.orig;
   cperm = col.perm;

   /*  move rhsidx to a heap
    */

   for(i = 0; i < rhsn;)
      enQueueMin(rhsidx, &i, cperm[rhsidx[i]]);

   ridx = u.row.idx;

   rval = u.row.val;

   rlen = u.row.len;

   rbeg = u.row.start;

   while(rhsn > 0)
   {
      i = deQueueMin(rhsidx, &rhsn);
      assert(i >= 0 && i < thedim);
      c = corig[i];
      assert(c >= 0 && c < thedim);
      x = rhs[c];
      rhs[c] = 0;

      if(isNotZero(x, eps))
      {
         r = rorig[i];
         assert(r >= 0 && r < thedim);
         x *= diag[r];
         vec[r] = x;
         k = rbeg[r];
         assert(k >= 0 && k < u.row.size);
         idx = &ridx[k];
         val = &rval[k];

         for(int m = rlen[r]; m; --m)
         {
            j = *idx++;
            assert(j >= 0 && j < thedim);
            y = rhs[j];

            if(y == 0)
            {
               y = -x * (*val++);

               if(isNotZero(y, eps))
               {
                  rhs[j] = y;
                  enQueueMin(rhsidx, &rhsn, cperm[j]);
               }
            }
            else
            {
               y -= x * (*val++);
               rhs[j] = (y != 0) ? y : SOPLEX_FACTOR_MARKER;
            }
         }
      }
   }
}


int CLUFactor::solveLleftForest(Real eps, Real* vec, int* nonz, int n)
{
   int i, j, k, end;
   Real x, y;
   Real* val, *lval;
   int* idx, *lidx, *lrow, *lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;
   end = l.firstUpdate;

   for(i = l.firstUnused - 1; i >= end; --i)
   {
      assert(i >= 0 && i < l.size);

      if((x = vec[lrow[i]]) != 0.0)
      {
         k = lbeg[i];
         assert(k >= 0 && k < l.size);
         val = &lval[k];
         idx = &lidx[k];

         for(j = lbeg[i + 1]; j > k; --j)
         {
            int m = *idx++;
            assert(m >= 0 && m < thedim);
            y = vec[m];

            if(y == 0)
            {
               y = -x * (*val++);

               if(isNotZero(y, eps))
               {
                  vec[m] = y;
                  nonz[n++] = m;
               }
            }
            else
            {
               y -= x * (*val++);
               vec[m] = (y != 0) ? y : SOPLEX_FACTOR_MARKER;
            }
         }
      }
   }

   return n;
}


void CLUFactor::solveLleftForestNoNZ(Real* vec)
{
   int i, j, k, end;
   Real x;
   Real* val, *lval;
   int* idx, *lidx, *lrow, *lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;
   end = l.firstUpdate;

   for(i = l.firstUnused - 1; i >= end; --i)
   {
      if((x = vec[lrow[i]]) != 0.0)
      {
         assert(i >= 0 && i < l.size);
         k = lbeg[i];
         assert(k >= 0 && k < l.size);
         val = &lval[k];
         idx = &lidx[k];

         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            vec[*idx++] -= x * (*val++);
         }
      }
   }
}


int CLUFactor::solveLleft(Real eps, Real* vec, int* nonz, int rn)
{
   int i, j, k, n;
   int r;
   Real x, y;
   Real* rval, *val;
   int* ridx, *idx;
   int* rbeg;
   int* rorig, *rperm;
   int* last;

   ridx  = l.ridx;
   rval  = l.rval;
   rbeg  = l.rbeg;
   rorig = l.rorig;
   rperm = l.rperm;
   n     = 0;

   i = l.firstUpdate - 1;
#ifndef WITH_L_ROWS
#pragma warn "Not yet implemented, define WITH_L_ROWS"
   Real*   lval = l.val;
   int*    lidx = l.idx;
   int*    lrow = l.row;
   int*    lbeg = l.start;

   for(; i >= 0; --i)
   {
      k   = lbeg[i];
      val = &lval[k];
      idx = &lidx[k];
      x   = 0;

      for(j = lbeg[i + 1]; j > k; --j)
         x += vec[*idx++] * (*val++);

      vec[lrow[i]] -= x;
   }

#else

   /*  move rhsidx to a heap
    */
   for(i = 0; i < rn;)
      enQueueMax(nonz, &i, rperm[nonz[i]]);

   last = nonz + thedim;

   while(rn > 0)
   {
      i = deQueueMax(nonz, &rn);
      r = rorig[i];
      x = vec[r];

      if(isNotZero(x, eps))
      {
         *(--last) = r;
         n++;
         k = rbeg[r];
         j = rbeg[r + 1] - k;
         val = &rval[k];
         idx = &ridx[k];

         while(j-- > 0)
         {
            assert(l.rperm[*idx] < i);
            int m = *idx++;
            y = vec[m];

            if(y == 0)
            {
               y = -x * *val++;

               if(isNotZero(y, eps))
               {
                  vec[m] = y;
                  enQueueMax(nonz, &rn, rperm[m]);
               }
            }
            else
            {
               y -= x * *val++;
               vec[m] = (y != 0) ? y : SOPLEX_FACTOR_MARKER;
            }
         }
      }
      else
         vec[r] = 0;
   }

   for(i = 0; i < n; ++i)
      *nonz++ = *last++;

#endif

   return n;
}


void CLUFactor::solveLleftNoNZ(Real* vec)
{
   int i, j, k;
   int r;
   Real x;
   Real* rval, *val;
   int* ridx, *idx;
   int* rbeg;
   int* rorig;

   ridx = l.ridx;
   rval = l.rval;
   rbeg = l.rbeg;
   rorig = l.rorig;

#ifndef WITH_L_ROWS
   Real* lval = l.val;
   int*    lidx = l.idx;
   int*    lrow = l.row;
   int*    lbeg = l.start;

   i = l.firstUpdate - 1;
   assert(i < thedim);

   for(; i >= 0; --i)
   {
      k = lbeg[i];
      assert(k >= 0 && k < l.size);
      val = &lval[k];
      idx = &lidx[k];
      x = 0;

      for(j = lbeg[i + 1]; j > k; --j)
      {
         assert(*idx >= 0 && *idx < thedim);
         x += vec[*idx++] * (*val++);
      }

      vec[lrow[i]] -= x;
   }

#else

   for(i = thedim; i--;)
   {
      r = rorig[i];
      x = vec[r];

      if(x != 0.0)
      {
         k = rbeg[r];
         j = rbeg[r + 1] - k;
         val = &rval[k];
         idx = &ridx[k];

         while(j-- > 0)
         {
            assert(l.rperm[*idx] < i);
            vec[*idx++] -= x * *val++;
         }
      }
   }

#endif
}

void inline CLUFactor::updateSolutionVectorLright(Real change, int j, Real& vec, int* idx, int& nnz)
{
   // create a new entry in #ridx
   if(vec == 0.0)
   {
      assert(nnz < thedim);
      idx[nnz] = j;
      ++nnz;
   }

   vec -= change;

   // mark the entry where exact eliminiation occurred
   if(vec == 0.0)
      vec = SOPLEX_FACTOR_MARKER;
}

// solve Lz = b, inplace, using and preserving sparisity structure in the rhs and solution vector
// arrays #vec and #ridx must be large enough to hold #thedim entries!
void CLUFactor::vSolveLright(Real* vec, int* ridx, int& rn, Real eps)
{
   int i, j, k, n;
   int end;
   Real x;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   // loop through columns of L
   for(i = 0; i < end; ++i)
   {
      x = vec[lrow[i]];

      // check whether there is a corresponding value in the rhs vector; skipping/ignoring FACTOR_MARKER
      if(isNotZero(x, eps))
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         // apply \f$- x * L_{k,i}\f$ to all corresponding values in rhs/solution vector
         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            n = *idx++;
            updateSolutionVectorLright(x * (*val), n, vec[n], ridx, rn);
            ++val;
         }
      }
   }

   if(l.updateType)                     /* Forest-Tomlin Updates */
   {
      end = l.firstUnused;

      for(; i < end; ++i)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         StableSum<Real> tmp;

         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            tmp += vec[*idx++] * (*val++);
         }

         j = lrow[i];
         x = tmp;

         if(isNotZero(x, eps))
            updateSolutionVectorLright(x, j, vec[j], ridx, rn);
      }
   }
}

// solve with L for two right hand sides
// see above methods for documentation
void CLUFactor::vSolveLright2(
   Real* vec, int* ridx, int& rn, Real eps,
   Real* vec2, int* ridx2, int& rn2, Real eps2)
{
   int i, j, k, n;
   int end;
   Real x, x2;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   // loop through columns of L
   for(i = 0; i < end; ++i)
   {
      j = lrow[i];
      x2 = vec2[j];
      x = vec[j];

      // check whether there is a corresponding value in the first rhs vector; skipping/ignoring FACTOR_MARKER
      if(isNotZero(x, eps))
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         // check whether there is  also a corresponding value in the second rhs vector; skipping/ignoring FACTOR_MARKER
         if(isNotZero(x2, eps2))
         {
            for(j = lbeg[i + 1]; j > k; --j)
            {
               assert(*idx >= 0 && *idx < thedim);
               n = *idx++;
               updateSolutionVectorLright(x * (*val), n, vec[n], ridx, rn);
               updateSolutionVectorLright(x2 * (*val), n, vec2[n], ridx2, rn2);
               ++val;
            }
         }
         // only the first vector needs to be modified
         else
         {
            for(j = lbeg[i + 1]; j > k; --j)
            {
               assert(*idx >= 0 && *idx < thedim);
               n = *idx++;
               updateSolutionVectorLright(x * (*val), n, vec[n], ridx, rn);
               ++val;
            }
         }
      }
      // only the second vector needs to be modified
      else if(isNotZero(x2, eps2))
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            n = *idx++;
            updateSolutionVectorLright(x2 * (*val), n, vec2[n], ridx2, rn2);
            ++val;
         }
      }
   }

   if(l.updateType)                     /* Forest-Tomlin Updates */
   {
      end = l.firstUnused;

      for(; i < end; ++i)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         StableSum<Real> tmp1;
         StableSum<Real> tmp2;

         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            tmp1 += vec[*idx] * (*val);
            tmp2 += vec2[*idx++] * (*val++);
         }

         x = tmp1;
         x2 = tmp2;

         j = lrow[i];

         if(isNotZero(x, eps))
            updateSolutionVectorLright(x, j, vec[j], ridx, rn);

         if(isNotZero(x2, eps2))
            updateSolutionVectorLright(x2, j, vec2[j], ridx2, rn2);
      }
   }
}

// solve with L for three right hand sides
// see above methods for documentation
void CLUFactor::vSolveLright3(
   Real* vec, int* ridx, int& rn, Real eps,
   Real* vec2, int* ridx2, int& rn2, Real eps2,
   Real* vec3, int* ridx3, int& rn3, Real eps3)
{
   int i, j, k, n;
   int end;
   Real x, x2, x3;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;

   end = l.firstUpdate;

   for(i = 0; i < end; ++i)
   {
      j = lrow[i];
      x = vec[j];
      x2 = vec2[j];
      x3 = vec3[j];

      if(isNotZero(x, eps))
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         if(isNotZero(x2, eps2))
         {
            if(isNotZero(x3, eps3))
            {
               // case 1: all three vectors are nonzero at j
               for(j = lbeg[i + 1]; j > k; --j)
               {
                  assert(*idx >= 0 && *idx < thedim);
                  n = *idx++;
                  updateSolutionVectorLright(x * (*val), n, vec[n], ridx, rn);
                  updateSolutionVectorLright(x2 * (*val), n, vec2[n], ridx2, rn2);
                  updateSolutionVectorLright(x3 * (*val), n, vec3[n], ridx3, rn3);
                  ++val;
               }
            }
            else
            {
               // case 2: 1 and 2 are nonzero at j
               for(j = lbeg[i + 1]; j > k; --j)
               {
                  assert(*idx >= 0 && *idx < thedim);
                  n = *idx++;
                  updateSolutionVectorLright(x * (*val), n, vec[n], ridx, rn);
                  updateSolutionVectorLright(x2 * (*val), n, vec2[n], ridx2, rn2);
                  ++val;
               }
            }
         }
         else if(isNotZero(x3, eps3))
         {
            // case 3: 1 and 3 are nonzero at j
            for(j = lbeg[i + 1]; j > k; --j)
            {
               assert(*idx >= 0 && *idx < thedim);
               n = *idx++;
               updateSolutionVectorLright(x * (*val), n, vec[n], ridx, rn);
               updateSolutionVectorLright(x3 * (*val), n, vec3[n], ridx3, rn3);
               ++val;
            }
         }
         else
         {
            // case 4: only 1 is nonzero at j
            for(j = lbeg[i + 1]; j > k; --j)
            {
               assert(*idx >= 0 && *idx < thedim);
               n = *idx++;
               updateSolutionVectorLright(x * (*val), n, vec[n], ridx, rn);
               ++val;
            }
         }
      }
      else if(isNotZero(x2, eps2))
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         if(isNotZero(x3, eps3))
         {
            // case 5: 2 and 3 are nonzero at j
            for(j = lbeg[i + 1]; j > k; --j)
            {
               assert(*idx >= 0 && *idx < thedim);
               n = *idx++;
               updateSolutionVectorLright(x2 * (*val), n, vec2[n], ridx2, rn2);
               updateSolutionVectorLright(x3 * (*val), n, vec3[n], ridx3, rn3);
               ++val;
            }
         }
         else
         {
            // case 6: only 2 is nonzero at j
            for(j = lbeg[i + 1]; j > k; --j)
            {
               assert(*idx >= 0 && *idx < thedim);
               n = *idx++;
               updateSolutionVectorLright(x2 * (*val), n, vec2[n], ridx2, rn2);
               ++val;
            }
         }
      }
      else if(isNotZero(x3, eps3))
      {
         // case 7: only 3 is nonzero at j
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            n = *idx++;
            updateSolutionVectorLright(x3 * (*val), n, vec3[n], ridx3, rn3);
            ++val;
         }
      }
   }

   if(l.updateType)                      /* Forest-Tomlin Updates */
   {
      end = l.firstUnused;

      for(; i < end; ++i)
      {
         k = lbeg[i];
         idx = &(lidx[k]);
         val = &(lval[k]);

         StableSum<Real> tmp1;
         StableSum<Real> tmp2;
         StableSum<Real> tmp3;

         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            tmp1 += vec[*idx] * (*val);
            tmp2 += vec2[*idx] * (*val);
            tmp3 += vec3[*idx++] * (*val++);
         }

         x = tmp1;
         x2 = tmp2;
         x3 = tmp3;

         j = lrow[i];

         if(isNotZero(x, eps))
            updateSolutionVectorLright(x, j, vec[j], ridx, rn);

         if(isNotZero(x2, eps2))
            updateSolutionVectorLright(x2, j, vec2[j], ridx2, rn2);

         if(isNotZero(x3, eps3))
            updateSolutionVectorLright(x3, j, vec3[j], ridx3, rn3);
      }
   }
}

int CLUFactor::vSolveUright(Real* vec, int* vidx,
                            Real* rhs, int* ridx, int rn, Real eps)
{
   int i, j, k, r, c, n;
   int* rorig, *corig;
   int* rperm;
   int* cidx, *clen, *cbeg;
   Real* cval;
   Real x, y;

   int* idx;
   Real* val;

   rorig = row.orig;
   corig = col.orig;
   rperm = row.perm;

   cidx = u.col.idx;
   cval = u.col.val;
   clen = u.col.len;
   cbeg = u.col.start;

   n = 0;

   while(rn > 0)
   {
      /*      Find nonzero with highest permuted row index and setup i and r
       */
      i = deQueueMax(ridx, &rn);
      assert(i >= 0 && i < thedim);
      r = rorig[i];
      assert(r >= 0 && r < thedim);

      x = diag[r] * rhs[r];
      rhs[r] = 0;

      if(isNotZero(x, eps))
      {
         c = corig[i];
         assert(c >= 0 && c < thedim);
         vidx[n++] = c;
         vec[c] = x;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
         {
            assert(*idx >= 0 && *idx < thedim);
            k = *idx++;
            assert(k >= 0 && k < thedim);
            y = rhs[k];

            if(y == 0)
            {
               y = -x * (*val++);

               if(isNotZero(y, eps))
               {
                  rhs[k] = y;
                  enQueueMax(ridx, &rn, rperm[k]);
               }
            }
            else
            {
               y -= x * (*val++);
               y += (y == 0) ? SOPLEX_FACTOR_MARKER : 0;
               rhs[k] = y;
            }
         }

         if(rn > i * verySparseFactor4right)
         {
            /* continue with dense case */
            for(i = *ridx; i >= 0; --i)
            {
               r = rorig[i];
               assert(r >= 0 && r < thedim);
               x = diag[r] * rhs[r];
               rhs[r] = 0;

               if(isNotZero(x, eps))
               {
                  c = corig[i];
                  assert(c >= 0 && c < thedim);
                  vidx[n++] = c;
                  vec[c] = x;
                  val = &cval[cbeg[c]];
                  idx = &cidx[cbeg[c]];
                  j = clen[c];

                  while(j-- > 0)
                  {
                     assert(*idx >= 0 && *idx < thedim);
                     rhs[*idx++] -= x * (*val++);
                  }
               }
            }

            break;
         }
      }
   }

   return n;
}


void CLUFactor::vSolveUrightNoNZ(Real* vec,
                                 Real* rhs, int* ridx, int rn, Real eps)
{
   int i, j, k, r, c;
   int* rorig, *corig;
   int* rperm;
   int* cidx, *clen, *cbeg;
   Real* cval;
   Real x, y;

   int* idx;
   Real* val;

   rorig = row.orig;
   corig = col.orig;
   rperm = row.perm;

   cidx = u.col.idx;
   cval = u.col.val;
   clen = u.col.len;
   cbeg = u.col.start;

   while(rn > 0)
   {
      if(rn > *ridx * verySparseFactor4right)
      {
         /* continue with dense case */
         for(i = *ridx; i >= 0; --i)
         {
            assert(i >= 0 && i < thedim);
            r = rorig[i];
            assert(r >= 0 && r < thedim);
            x = diag[r] * rhs[r];
            rhs[r] = 0;

            if(isNotZero(x, eps))
            {
               c = corig[i];
               vec[c] = x;
               val = &cval[cbeg[c]];
               idx = &cidx[cbeg[c]];
               j = clen[c];

               while(j-- > 0)
               {
                  assert(*idx >= 0 && *idx < thedim);
                  rhs[*idx++] -= x * (*val++);
               }
            }
         }

         break;
      }

      /*      Find nonzero with highest permuted row index and setup i and r
       */
      i = deQueueMax(ridx, &rn);

      assert(i >= 0 && i < thedim);

      r = rorig[i];

      assert(r >= 0 && r < thedim);

      x = diag[r] * rhs[r];

      rhs[r] = 0;

      if(isNotZero(x, eps))
      {
         c = corig[i];
         vec[c] = x;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
         {
            k = *idx++;
            assert(k >= 0 && k < thedim);
            y = rhs[k];

            if(y == 0)
            {
               y = -x * (*val++);

               if(isNotZero(y, eps))
               {
                  rhs[k] = y;
                  enQueueMax(ridx, &rn, rperm[k]);
               }
            }
            else
            {
               y -= x * (*val++);
               y += (y == 0) ? SOPLEX_FACTOR_MARKER : 0;
               rhs[k] = y;
            }
         }
      }
   }
}


int CLUFactor::vSolveUright2(
   Real* vec, int* vidx, Real* rhs, int* ridx, int rn, Real eps,
   Real* vec2, Real* rhs2, int* ridx2, int rn2, Real eps2)
{
   int i, j, k, r, c, n;
   int* rorig, *corig;
   int* rperm;
   int* cidx, *clen, *cbeg;
   Real* cval;
   Real x, y;
   Real x2, y2;

   int* idx;
   Real* val;

   rorig = row.orig;
   corig = col.orig;
   rperm = row.perm;

   cidx = u.col.idx;
   cval = u.col.val;
   clen = u.col.len;
   cbeg = u.col.start;

   n = 0;

   while(rn + rn2 > 0)
   {
      /*      Find nonzero with highest permuted row index and setup i and r
       */
      if(rn <= 0)
         i = deQueueMax(ridx2, &rn2);
      else if(rn2 <= 0)
         i = deQueueMax(ridx, &rn);
      else if(*ridx2 > *ridx)
         i = deQueueMax(ridx2, &rn2);
      else if(*ridx2 < *ridx)
         i = deQueueMax(ridx, &rn);
      else
      {
         (void) deQueueMax(ridx, &rn);
         i = deQueueMax(ridx2, &rn2);
      }

      assert(i >= 0 && i < thedim);

      r = rorig[i];
      assert(r >= 0 && r < thedim);

      x = diag[r] * rhs[r];
      x2 = diag[r] * rhs2[r];
      rhs[r] = 0;
      rhs2[r] = 0;

      if(isNotZero(x, eps))
      {
         c = corig[i];
         vidx[n++] = c;
         vec[c] = x;
         vec2[c] = x2;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         if(isNotZero(x2, eps2))
         {
            while(j-- > 0)
            {
               k = *idx++;
               assert(k >= 0 && k < thedim);
               y2 = rhs2[k];

               if(y2 == 0)
               {
                  y2 = -x2 * (*val);

                  if(isNotZero(y2, eps2))
                  {
                     rhs2[k] = y2;
                     enQueueMax(ridx2, &rn2, rperm[k]);
                  }
               }
               else
               {
                  y2 -= x2 * (*val);
                  rhs2[k] = (y2 != 0) ? y2 : SOPLEX_FACTOR_MARKER;
               }

               y = rhs[k];

               if(y == 0)
               {
                  y = -x * (*val++);

                  if(isNotZero(y, eps))
                  {
                     rhs[k] = y;
                     enQueueMax(ridx, &rn, rperm[k]);
                  }
               }
               else
               {
                  y -= x * (*val++);
                  y += (y == 0) ? SOPLEX_FACTOR_MARKER : 0;
                  rhs[k] = y;
               }
            }
         }
         else
         {
            while(j-- > 0)
            {
               k = *idx++;
               assert(k >= 0 && k < thedim);
               y = rhs[k];

               if(y == 0)
               {
                  y = -x * (*val++);

                  if(isNotZero(y, eps))
                  {
                     rhs[k] = y;
                     enQueueMax(ridx, &rn, rperm[k]);
                  }
               }
               else
               {
                  y -= x * (*val++);
                  y += (y == 0) ? SOPLEX_FACTOR_MARKER : 0;
                  rhs[k] = y;
               }
            }
         }
      }
      else if(isNotZero(x2, eps2))
      {
         c = corig[i];
         assert(c >= 0 && c < thedim);
         vec2[c] = x2;
         val = &cval[cbeg[c]];
         idx = &cidx[cbeg[c]];
         j = clen[c];

         while(j-- > 0)
         {
            k = *idx++;
            assert(k >= 0 && k < thedim);
            y2 = rhs2[k];

            if(y2 == 0)
            {
               y2 = -x2 * (*val++);

               if(isNotZero(y2, eps2))
               {
                  rhs2[k] = y2;
                  enQueueMax(ridx2, &rn2, rperm[k]);
               }
            }
            else
            {
               y2 -= x2 * (*val++);
               rhs2[k] = (y2 != 0) ? y2 : SOPLEX_FACTOR_MARKER;
            }
         }
      }

      if(rn + rn2 > i * verySparseFactor4right)
      {
         /* continue with dense case */
         if(*ridx > *ridx2)
            i = *ridx;
         else
            i = *ridx2;

         for(; i >= 0; --i)
         {
            assert(i < thedim);
            r = rorig[i];
            assert(r >= 0 && r < thedim);
            x = diag[r] * rhs[r];
            x2 = diag[r] * rhs2[r];
            rhs[r] = 0;
            rhs2[r] = 0;

            if(isNotZero(x2, eps2))
            {
               c = corig[i];
               assert(c >= 0 && c < thedim);
               vec2[c] = x2;
               val = &cval[cbeg[c]];
               idx = &cidx[cbeg[c]];
               j = clen[c];

               if(isNotZero(x, eps))
               {
                  vidx[n++] = c;
                  vec[c] = x;

                  while(j-- > 0)
                  {
                     assert(*idx >= 0 && *idx < thedim);
                     rhs[*idx] -= x * (*val);
                     rhs2[*idx++] -= x2 * (*val++);
                  }
               }
               else
               {
                  while(j-- > 0)
                  {
                     assert(*idx >= 0 && *idx < thedim);
                     rhs2[*idx++] -= x2 * (*val++);
                  }
               }
            }
            else if(isNotZero(x, eps))
            {
               c = corig[i];
               assert(c >= 0 && c < thedim);
               vidx[n++] = c;
               vec[c] = x;
               val = &cval[cbeg[c]];
               idx = &cidx[cbeg[c]];
               j = clen[c];

               while(j-- > 0)
               {
                  assert(*idx >= 0 && *idx < thedim);
                  rhs[*idx++] -= x * (*val++);
               }
            }
         }

         break;
      }
   }

   return n;
}

int CLUFactor::vSolveUpdateRight(Real* vec, int* ridx, int n, Real eps)
{
   int i, j, k;
   int end;
   Real x, y;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   assert(!l.updateType);               /* no Forest-Tomlin Updates */

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;
   end = l.firstUnused;

   for(i = l.firstUpdate; i < end; ++i)
   {
      assert(i >= 0 && i < thedim);
      x = vec[lrow[i]];

      if(isNotZero(x, eps))
      {
         k = lbeg[i];
         assert(k >= 0 && k < l.size);
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            int m = ridx[n] = *idx++;
            assert(m >= 0 && m < thedim);
            y = vec[m];
            n += (y == 0) ? 1 : 0;
            y = y - x * (*val++);
            vec[m] = (y != 0) ? y : SOPLEX_FACTOR_MARKER;
         }
      }
   }

   return n;
}

void CLUFactor::vSolveUpdateRightNoNZ(Real* vec, Real /*eps*/)
{
   int i, j, k;
   int end;
   Real x;
   Real* lval, *val;
   int* lrow, *lidx, *idx;
   int* lbeg;

   assert(!l.updateType);               /* no Forest-Tomlin Updates */

   lval = l.val;
   lidx = l.idx;
   lrow = l.row;
   lbeg = l.start;
   end = l.firstUnused;

   for(i = l.firstUpdate; i < end; ++i)
   {
      assert(i >= 0 && i < thedim);

      if((x = vec[lrow[i]]) != 0.0)
      {
         k = lbeg[i];
         assert(k >= 0 && k < l.size);
         idx = &(lidx[k]);
         val = &(lval[k]);

         for(j = lbeg[i + 1]; j > k; --j)
         {
            assert(*idx >= 0 && *idx < thedim);
            vec[*idx++] -= x * (*val++);
         }
      }
   }
}


int CLUFactor::vSolveRight4update(Real eps,
                                  Real* vec, int* idx,                       /* result */
                                  Real* rhs, int* ridx, int rn,              /* rhs    */
                                  Real* forest, int* forestNum, int* forestIdx)
{
   vSolveLright(rhs, ridx, rn, eps);
   assert(rn >= 0 && rn <= thedim);

   /*  turn index list into a heap
    */

   if(forest)
   {
      Real x;
      int i, j, k;
      int* rperm;
      int* it = forestIdx;

      rperm = row.perm;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
         {
            enQueueMax(ridx, &j, rperm[*it++ = k]);
            forest[k] = x;
         }
         else
            rhs[k] = 0;
      }

      *forestNum = rn = j;
   }
   else
   {
      Real x;
      int i, j, k;
      int* rperm;

      rperm = row.perm;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
            enQueueMax(ridx, &j, rperm[k]);
         else
            rhs[k] = 0;
      }

      rn = j;
   }

   rn = vSolveUright(vec, idx, rhs, ridx, rn, eps);

   if(!l.updateType)             /* no Forest-Tomlin Updates */
      rn = vSolveUpdateRight(vec, idx, rn, eps);

   return rn;
}

int CLUFactor::vSolveRight4update2(Real eps,
                                   Real* vec, int* idx,                  /* result1 */
                                   Real* rhs, int* ridx, int rn,         /* rhs1    */
                                   Real* vec2, Real eps2,              /* result2 */
                                   Real* rhs2, int* ridx2, int rn2,      /* rhs2    */
                                   Real* forest, int* forestNum, int* forestIdx)
{
   vSolveLright2(rhs, ridx, rn, eps, rhs2, ridx2, rn2, eps2);
   assert(rn >= 0 && rn <= thedim);
   assert(rn2 >= 0 && rn2 <= thedim);

   /*  turn index list into a heap
    */

   if(forest)
   {
      Real x;
      int i, j, k;
      int* rperm;
      int* it = forestIdx;

      rperm = row.perm;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
         {
            enQueueMax(ridx, &j, rperm[*it++ = k]);
            forest[k] = x;
         }
         else
            rhs[k] = 0;
      }

      *forestNum = rn = j;
   }
   else
   {
      Real x;
      int i, j, k;
      int* rperm;

      rperm = row.perm;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
            enQueueMax(ridx, &j, rperm[k]);
         else
            rhs[k] = 0;
      }

      rn = j;
   }

   if(rn2 > thedim * verySparseFactor4right)
   {
      ridx2[0] = thedim - 1;
      /* ridx2[1] = thedim - 2; */
   }
   else
   {
      Real x;
      /*      Real  maxabs; */
      int i, j, k;
      int* rperm;

      /*      maxabs = 1;    */
      rperm = row.perm;

      for(i = j = 0; i < rn2; ++i)
      {
         k = ridx2[i];
         assert(k >= 0 && k < thedim);
         x = rhs2[k];

         if(x < -eps2)
         {
            /*              maxabs = (maxabs < -x) ? -x : maxabs;  */
            enQueueMax(ridx2, &j, rperm[k]);
         }
         else if(x > eps2)
         {
            /*              maxabs = (maxabs < x) ? x : maxabs;    */
            enQueueMax(ridx2, &j, rperm[k]);
         }
         else
            rhs2[k] = 0;
      }

      rn2 = j;

      /*      eps2 = maxabs * eps2;  */
   }

   rn = vSolveUright(vec, idx, rhs, ridx, rn, eps);

   vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2, eps2);

   /*
    *  rn = vSolveUright2(vec, idx, rhs, ridx, rn, eps, vec2, rhs2, ridx2, rn2, eps2);
    */

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = vSolveUpdateRight(vec, idx, rn, eps);
      vSolveUpdateRightNoNZ(vec2, eps2);
   }

   return rn;
}

void CLUFactor::vSolveRight4update2sparse(Real eps, Real* vec, int* idx,         /* result1 */
      Real* rhs, int* ridx, int& rn,        /* rhs1    */
      Real eps2, Real* vec2, int* idx2,     /* result2 */
      Real* rhs2, int* ridx2, int& rn2,     /* rhs2    */
      Real* forest, int* forestNum, int* forestIdx)
{
   /* solve with L */
   vSolveLright2(rhs, ridx, rn, eps, rhs2, ridx2, rn2, eps2);
   assert(rn >= 0 && rn <= thedim);
   assert(rn2 >= 0 && rn2 <= thedim);

   Real x;
   int i, j, k;
   int* rperm = row.perm;

   /*  turn index list into a heap for both ridx and ridx2 */
   if(forest)
   {
      int* it = forestIdx;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
         {
            enQueueMax(ridx, &j, rperm[*it++ = k]);
            forest[k] = x;
         }
         else
            rhs[k] = 0;
      }

      *forestNum = rn = j;
   }
   else
   {
      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
            enQueueMax(ridx, &j, rperm[k]);
         else
            rhs[k] = 0;
      }

      rn = j;
   }

   for(i = j = 0; i < rn2; ++i)
   {
      k = ridx2[i];
      assert(k >= 0 && k < thedim);
      x = rhs2[k];

      if(isNotZero(x, eps2))
         enQueueMax(ridx2, &j, rperm[k]);
      else
         rhs2[k] = 0;
   }

   rn2 = j;

   /* solve with U */
   rn = vSolveUright(vec, idx, rhs, ridx, rn, eps);
   rn2 = vSolveUright(vec2, idx2, rhs2, ridx2, rn2, eps2);

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = vSolveUpdateRight(vec, idx, rn, eps);
      rn2 = vSolveUpdateRight(vec2, idx2, rn2, eps2);
   }
}


int CLUFactor::vSolveRight4update3(Real eps,
                                   Real* vec, int* idx,                 /* result1 */
                                   Real* rhs, int* ridx, int rn,        /* rhs1    */
                                   Real* vec2, Real eps2,               /* result2 */
                                   Real* rhs2, int* ridx2, int rn2,     /* rhs2    */
                                   Real* vec3, Real eps3,               /* result3 */
                                   Real* rhs3, int* ridx3, int rn3,     /* rhs3    */
                                   Real* forest, int* forestNum, int* forestIdx)
{

   vSolveLright3(rhs, ridx, rn, eps, rhs2, ridx2, rn2, eps2, rhs3, ridx3, rn3, eps3);
   assert(rn >= 0 && rn <= thedim);
   assert(rn2 >= 0 && rn2 <= thedim);
   assert(rn3 >= 0 && rn3 <= thedim);

   /*  turn index list into a heap
    */

   if(forest)
   {
      Real x;
      int i, j, k;
      int* rperm;
      int* it = forestIdx;

      rperm = row.perm;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
         {
            enQueueMax(ridx, &j, rperm[*it++ = k]);
            forest[k] = x;
         }
         else
            rhs[k] = 0;
      }

      *forestNum = rn = j;
   }
   else
   {
      Real x;
      int i, j, k;
      int* rperm;

      rperm = row.perm;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
            enQueueMax(ridx, &j, rperm[k]);
         else
            rhs[k] = 0;
      }

      rn = j;
   }

   if(rn2 > thedim * verySparseFactor4right)
   {
      ridx2[0] = thedim - 1;
   }
   else
   {
      Real x;
      int i, j, k;
      int* rperm;

      rperm = row.perm;

      for(i = j = 0; i < rn2; ++i)
      {
         k = ridx2[i];
         assert(k >= 0 && k < thedim);
         x = rhs2[k];

         if(x < -eps2)
         {
            enQueueMax(ridx2, &j, rperm[k]);
         }
         else if(x > eps2)
         {
            enQueueMax(ridx2, &j, rperm[k]);
         }
         else
            rhs2[k] = 0;
      }

      rn2 = j;
   }

   if(rn3 > thedim * verySparseFactor4right)
   {
      ridx3[0] = thedim - 1;
   }
   else
   {
      Real x;
      int i, j, k;
      int* rperm;

      rperm = row.perm;

      for(i = j = 0; i < rn3; ++i)
      {
         k = ridx3[i];
         assert(k >= 0 && k < thedim);
         x = rhs3[k];

         if(x < -eps3)
         {
            enQueueMax(ridx3, &j, rperm[k]);
         }
         else if(x > eps3)
         {
            enQueueMax(ridx3, &j, rperm[k]);
         }
         else
            rhs3[k] = 0;
      }

      rn3 = j;
   }

   rn = vSolveUright(vec, idx, rhs, ridx, rn, eps);

   vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2, eps2);
   vSolveUrightNoNZ(vec3, rhs3, ridx3, rn3, eps3);

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = vSolveUpdateRight(vec, idx, rn, eps);
      vSolveUpdateRightNoNZ(vec2, eps2);
      vSolveUpdateRightNoNZ(vec3, eps3);
   }

   return rn;
}

void CLUFactor::vSolveRight4update3sparse(Real eps, Real* vec, int* idx,         /* result1 */
      Real* rhs, int* ridx, int& rn,        /* rhs1    */
      Real eps2, Real* vec2, int* idx2,     /* result2 */
      Real* rhs2, int* ridx2, int& rn2,     /* rhs2    */
      Real eps3, Real* vec3, int* idx3,     /* result3 */
      Real* rhs3, int* ridx3, int& rn3,     /* rhs3    */
      Real* forest, int* forestNum, int* forestIdx)
{
   vSolveLright3(rhs, ridx, rn, eps, rhs2, ridx2, rn2, eps2, rhs3, ridx3, rn3, eps3);
   assert(rn >= 0 && rn <= thedim);
   assert(rn2 >= 0 && rn2 <= thedim);
   assert(rn3 >= 0 && rn3 <= thedim);

   Real x;
   int i, j, k;
   int* rperm = row.perm;

   /*  turn index list into a heap */
   if(forest)
   {
      int* it = forestIdx;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
         {
            enQueueMax(ridx, &j, rperm[*it++ = k]);
            forest[k] = x;
         }
         else
            rhs[k] = 0;
      }

      *forestNum = rn = j;
   }
   else
   {
      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(isNotZero(x, eps))
            enQueueMax(ridx, &j, rperm[k]);
         else
            rhs[k] = 0;
      }

      rn = j;
   }

   for(i = j = 0; i < rn2; ++i)
   {
      k = ridx2[i];
      assert(k >= 0 && k < thedim);
      x = rhs2[k];

      if(isNotZero(x, eps2))
         enQueueMax(ridx2, &j, rperm[k]);
      else
         rhs2[k] = 0;
   }

   rn2 = j;

   for(i = j = 0; i < rn3; ++i)
   {
      k = ridx3[i];
      assert(k >= 0 && k < thedim);
      x = rhs3[k];

      if(isNotZero(x, eps3))
         enQueueMax(ridx3, &j, rperm[k]);
      else
         rhs3[k] = 0;
   }

   rn3 = j;

   rn = vSolveUright(vec, idx, rhs, ridx, rn, eps);
   rn2 = vSolveUright(vec2, idx2, rhs2, ridx2, rn2, eps2);
   rn3 = vSolveUright(vec3, idx3, rhs3, ridx3, rn3, eps3);

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = vSolveUpdateRight(vec, idx, rn, eps);
      rn2 = vSolveUpdateRight(vec2, idx2, rn2, eps2);
      rn3 = vSolveUpdateRight(vec3, idx3, rn3, eps3);
   }
}

void CLUFactor::vSolveRightNoNZ(
   Real* vec, Real eps,             /* result */
   Real* rhs, int* ridx, int rn)    /* rhs    */
{
   vSolveLright(rhs, ridx, rn, eps);
   assert(rn >= 0 && rn <= thedim);

   if(rn > thedim * verySparseFactor4right)
   {
      *ridx = thedim - 1;
   }
   else
   {
      Real x;
      /*      Real  maxabs; */
      int i, j, k;
      int* rperm;

      /*      maxabs = 1;    */
      rperm = row.perm;

      for(i = j = 0; i < rn; ++i)
      {
         k = ridx[i];
         assert(k >= 0 && k < thedim);
         x = rhs[k];

         if(x < -eps)
         {
            /*              maxabs = (maxabs < -x) ? -x : maxabs;  */
            enQueueMax(ridx, &j, rperm[k]);
         }
         else if(x > eps)
         {
            /*              maxabs = (maxabs < x) ? x : maxabs;    */
            enQueueMax(ridx, &j, rperm[k]);
         }
         else
            rhs[k] = 0;
      }

      rn = j;

      /*      eps2 = maxabs * eps2;  */
   }

   vSolveUrightNoNZ(vec, rhs, ridx, rn, eps);

   if(!l.updateType)             /* no Forest-Tomlin Updates */
      vSolveUpdateRightNoNZ(vec, eps);
}

int CLUFactor::vSolveLeft(Real eps,
                          Real* vec, int* idx,                       /* result */
                          Real* rhs, int* ridx, int rn)            /* rhs    */
{

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = solveUpdateLeft(eps, rhs, ridx, rn);
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
   }
   else
   {
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn = solveLleftForest(eps, vec, idx, rn);
   }

   // TODO verify the correctness of this check
   if(rn + l.firstUpdate > verySparseFactor4left * thedim)
   {
      // perform the dense solve
      solveLleftNoNZ(vec);
      // signal the caller that the nonzero pattern is lost
      return 0;
   }
   else
      return solveLleft(eps, vec, idx, rn);
}

int CLUFactor::vSolveLeft2(Real eps,
                           Real* vec, int* idx,                      /* result */
                           Real* rhs, int* ridx, int rn,             /* rhs    */
                           Real* vec2,                               /* result2 */
                           Real* rhs2, int* ridx2, int rn2)        /* rhs2    */
{

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = solveUpdateLeft(eps, rhs, ridx, rn);
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn2 = solveUpdateLeft(eps, rhs2, ridx2, rn2);
      solveUleftNoNZ(eps, vec2, rhs2, ridx2, rn2);
   }
   else
   {
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn = solveLleftForest(eps, vec, idx, rn);
      solveUleftNoNZ(eps, vec2, rhs2, ridx2, rn2);
      solveLleftForestNoNZ(vec2);
   }

   rn = solveLleft(eps, vec, idx, rn);

   solveLleftNoNZ(vec2);

   return rn;
}

void CLUFactor::vSolveLeft2sparse(Real eps,
                                  Real* vec, int* idx,                      /* result */
                                  Real* rhs, int* ridx, int& rn,            /* rhs    */
                                  Real* vec2, int* idx2,                    /* result2 */
                                  Real* rhs2, int* ridx2, int& rn2)         /* rhs2    */
{
   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = solveUpdateLeft(eps, rhs, ridx, rn);
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn2 = solveUpdateLeft(eps, rhs2, ridx2, rn2);
      rn2 = solveUleft(eps, vec2, idx2, rhs2, ridx2, rn2);
   }
   else
   {
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn = solveLleftForest(eps, vec, idx, rn);
      rn2 = solveUleft(eps, vec2, idx2, rhs2, ridx2, rn2);
      rn2 = solveLleftForest(eps, vec2, idx2, rn2);

   }

   rn = solveLleft(eps, vec, idx, rn);
   rn2 = solveLleft(eps, vec2, idx2, rn2);
}


int CLUFactor::vSolveLeft3(Real eps,
                           Real* vec, int* idx,                      /* result */
                           Real* rhs, int* ridx, int rn,             /* rhs    */
                           Real* vec2,                               /* result2 */
                           Real* rhs2, int* ridx2, int rn2,          /* rhs2    */
                           Real* vec3,                               /* result3 */
                           Real* rhs3, int* ridx3, int rn3)          /* rhs3    */
{

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = solveUpdateLeft(eps, rhs, ridx, rn);
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn2 = solveUpdateLeft(eps, rhs2, ridx2, rn2);
      solveUleftNoNZ(eps, vec2, rhs2, ridx2, rn2);
      rn3 = solveUpdateLeft(eps, rhs3, ridx3, rn3);
      solveUleftNoNZ(eps, vec3, rhs3, ridx3, rn3);
   }
   else
   {
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn = solveLleftForest(eps, vec, idx, rn);
      solveUleftNoNZ(eps, vec2, rhs2, ridx2, rn2);
      solveLleftForestNoNZ(vec2);
      solveUleftNoNZ(eps, vec3, rhs3, ridx3, rn3);
      solveLleftForestNoNZ(vec3);
   }

   rn = solveLleft(eps, vec, idx, rn);

   solveLleftNoNZ(vec2);
   solveLleftNoNZ(vec3);

   return rn;
}

void CLUFactor::vSolveLeft3sparse(Real eps,
                                  Real* vec, int* idx,                      /* result */
                                  Real* rhs, int* ridx, int& rn,            /* rhs    */
                                  Real* vec2, int* idx2,                    /* result2 */
                                  Real* rhs2, int* ridx2, int& rn2,         /* rhs2    */
                                  Real* vec3, int* idx3,                    /* result3 */
                                  Real* rhs3, int* ridx3, int& rn3)         /* rhs3    */
{
   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn = solveUpdateLeft(eps, rhs, ridx, rn);
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn2 = solveUpdateLeft(eps, rhs2, ridx2, rn2);
      rn2 = solveUleft(eps, vec2, idx2, rhs2, ridx2, rn2);
      rn3 = solveUpdateLeft(eps, rhs3, ridx3, rn3);
      rn3 = solveUleft(eps, vec3, idx3, rhs3, ridx3, rn3);
   }
   else
   {
      rn = solveUleft(eps, vec, idx, rhs, ridx, rn);
      rn = solveLleftForest(eps, vec, idx, rn);
      rn2 = solveUleft(eps, vec2, idx2, rhs2, ridx2, rn2);
      rn2 = solveLleftForest(eps, vec2, idx2, rn2);
      rn3 = solveUleft(eps, vec3, idx3, rhs3, ridx3, rn3);
      rn3 = solveLleftForest(eps, vec3, idx3, rn3);
   }

   rn = solveLleft(eps, vec, idx, rn);
   rn2 = solveLleft(eps, vec2, idx2, rn2);
   rn3 = solveLleft(eps, vec3, idx3, rn3);
}


void CLUFactor::vSolveLeftNoNZ(Real eps,
                               Real* vec2,                            /* result2 */
                               Real* rhs2, int* ridx2, int rn2)     /* rhs2    */
{

   if(!l.updateType)             /* no Forest-Tomlin Updates */
   {
      rn2 = solveUpdateLeft(eps, rhs2, ridx2, rn2);
      solveUleftNoNZ(eps, vec2, rhs2, ridx2, rn2);
   }
   else
   {
      solveUleftNoNZ(eps, vec2, rhs2, ridx2, rn2);
      solveLleftForestNoNZ(vec2);
   }

   solveLleftNoNZ(vec2);
}
} // namespace soplex