450-soplex Report

**450-soplex**
Exec Time (%)	Codelet Name	Error (%)
47.4	spxsteeppr__ZN6soplex10SPxSteepPR9enter...	98.1
12.0	spxsteeppr__ZN6soplex10SPxSteepPR12sele...	5.2
4.7	ssvector__ZN6soplex8SSVector5setupEv_212	5.2
8.1	enter__ZN6soplex6SoPlex10updateTestEv_2...	1.4
1.9	ssvector__ZN6soplex8SSVector5clearEv_74	3.0
15.1	ssvector__ZN6soplex8SSVector18assign2pr...	5.3

Exec Time : 8.433491e+11 cycles

Matching : 39.9%

                /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1997-1999 Roland Wunderling                              */
/*                  1997-2002 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.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef SPEC_CPU
#pragma ident "@(#) $Id: spxsteeppr.cpp,v 1.18 2002/03/21 16:06:19 bzfkocht Exp $"
#endif

//#define DEBUGGING 1

#include <assert.h>
#include <iostream>

#include "spxdefines.h"
#include "message.h"
#include "spxsteeppr.h"
#include "random.h"

namespace soplex
{

#define EQ_PREF 1000

void SPxSteepPR::clear()
{
   thesolver = 0;
   prefSetup = 0;
}

void SPxSteepPR::load(SoPlex* base)
{
   thesolver = base;

   if (base)
   {
      workVec.clear();
      workVec.reDim(base->dim());
      workRhs.clear();
      workRhs.reDim(base->dim());

      leavePref.reSize(base->dim());
      coPref.reSize (base->dim());
      pref.reSize (base->coDim());
      prefSetup = 0;
   }
}

void SPxSteepPR::setType(SoPlex::Type type)
{
   int i;

   workRhs.epsilon = thesolver->epsilon();

   pref.reSize (thesolver->coDim());
   coPref.reSize(thesolver->dim());
   setupPrefs(type);

   if (setup == DEFAULT)
   {
      if (type == SoPlex::ENTER)
      {
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
            // coPenalty[i] = 10;
            coPenalty[i] = 2;
         penalty.reDim(thesolver->coDim());
         for (i = thesolver->coDim() - 1; i >= 0; --i)
            // penalty[i] = 10;
            penalty[i] = 1;
         // penalty[i] = 1 + thesolver->vector(i).size() / thesolver->dim();
      }
      else
      {
         assert(type == SoPlex::LEAVE);
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
         {
            // coPenalty[i] = 1;
            SPxId id = thesolver->basis().baseId(i);
            int n = thesolver->number(id);
            if (thesolver->isId(id))
               leavePref[i] = pref[n];
            else
               leavePref[i] = coPref[n];
            coPenalty[i] = 1 + thesolver->basis().baseVec(i).size()
                           / Real(thesolver->dim());
         }
      }
   }
   else
   {
      std::cerr << 
         "sorry, no exact setup for steepest edge multipliers implemented\n";

      if (type == SoPlex::ENTER)
      {
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
            coPenalty[i] = 1;
         penalty.reDim(thesolver->coDim());
         for (i = thesolver->coDim() - 1; i >= 0; --i)
            penalty[i] = 1 + thesolver->vector(i).length2();
      }
      else
      {
         assert(type == SoPlex::LEAVE);
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
         {
            coPenalty[i] = 1 + thesolver->basis().baseVec(i).size()
                           / Real(thesolver->dim());
         }
      }
   }

   workVec.clear();
   workRhs.clear();
}

void SPxSteepPR::setupPrefs(Real mult, Real /*tie*/, Real /*cotie*/,
                             Real shift, Real coshift,
                             int rs, int cs, int re, int ce)
{
   Real *p;
   Real *cp;
   Real *end;
   // Real rtie, ctie;
   Real rshift, cshift;
   int i;

   if (thesolver->rep() == SoPlex::COLUMN)
   {
      cp = pref.get_ptr();
      p = coPref.get_ptr();
      // ctie = tie;
      // rtie = cotie;
      cshift = shift;
      rshift = coshift;
   }
   else
   {
      p = pref.get_ptr();
      cp = coPref.get_ptr();
      // rtie = tie;
      // ctie = cotie;
      rshift = shift;
      cshift = coshift;
   }

   if (re < 0)
      re = thesolver->nRows();
   for (i = re; --i >= rs;)
   {
      p[i] = rshift;
      //      p[i] += rtie * thesolver->rowVector(i).size() / Real(thesolver->nCols());
      //      p[i] += EQ_PREF * (thesolver->rhs(i) == thesolver->lhs(i));
      //      p[i] += EQ_PREF * (thesolver->rhs(i) >=  infinity
      //                     &&  thesolver->lhs(i) <= -infinity);
   }

   if (ce < 0)
      ce = thesolver->nCols();
   for (i = ce; --i >= cs;)
   {
      cp[i] = cshift;
      //      cp[i] += ctie * thesolver->colVector(i).size() / Real(thesolver->nRows());
      //      cp[i] += EQ_PREF * (thesolver->upper(i) == thesolver->lower(i));
      //      cp[i] += EQ_PREF * (thesolver->upper(i) >=  infinity
      //                      &&  thesolver->lower(i) <= -infinity);
   }


   i = 0;
   cp = coPref.get_ptr();
   end = cp + coPref.size();
   while (cp < end)
      *cp++ *= 1.0 - mult * i++;

   p = pref.get_ptr();
   end = p + pref.size();
   i = pref.size();
   while (p < end)
      *p++ *= 1.0 + mult * i--;
}

void SPxSteepPR::setupPrefs(SoPlex::Type tp)
{
   if (tp != prefSetup)
   {
      Real mult = 1e-8 / Real(1 + thesolver->dim() + thesolver->coDim());
      if (tp == SoPlex::ENTER)
         setupPrefs(-mult, -1e-5, -1e-5, 1, 1);
      else
         setupPrefs(mult, 1e-5, 1e-5, 1, 1);
      prefSetup = tp;
   }
}

void SPxSteepPR::setRep(SoPlex::Representation)
{
   if (workVec.dim() != thesolver->dim())
   {
      DVector tmp = penalty;
      penalty = coPenalty;
      coPenalty = tmp;

      workVec.clear();
      workVec.reDim(thesolver->dim());
   }
}

void SPxSteepPR::left4X(int n, SPxId id, int start, int incr)
{
   assert(thesolver->type() == SoPlex::LEAVE);

   if (id.isValid())
   {
      int i, j;
      Real x;
      // Real               delta         = 0.1;   // thesolver->epsilon();
      Real delta = 0.1 + 1.0 / thesolver->basis().iteration();
      Real* coPenalty_ptr = coPenalty.get_ptr();
      const Real* workVec_ptr = workVec.get_const_ptr();
      const Real* rhoVec = thesolver->fVec().delta().values();
      Real rhov_1 = 1 / rhoVec[n];
      Real beta_q = thesolver->coPvec().delta().length2()
                      * rhov_1 * rhov_1;

      assert(rhoVec[n] >= theeps || -rhoVec[n] >= theeps);

      //  Update #coPenalty# vector
      const IdxSet& rhoIdx = thesolver->fVec().idx();
      int len = thesolver->fVec().idx().size();
      for (i = len - 1 - start; i >= 0; i -= incr)
      {
         j = rhoIdx.index(i);
         x = coPenalty_ptr[j] += rhoVec[j]
            * (beta_q * rhoVec[j] - 2 * rhov_1 * workVec_ptr[j]);
         if (x < delta)
            // coPenalty_ptr[j] = delta / (1+delta-x);
            coPenalty_ptr[j] = delta;
         else if (x >= infinity)
            coPenalty_ptr[j] = 1 / theeps;
      }

      coPenalty_ptr[n] = beta_q;
      //@ coPenalty_ptr[n] = 0.999*beta_q;
      //@ coPenalty_ptr[n] = 1.001*beta_q;
   }

}

void SPxSteepPR::left4(int n, SPxId id)
{
   //  Update preference multiplier in #leavePref#
   if (thesolver->isId(id))
      leavePref[n] = pref[thesolver->number(id)];
   else if (thesolver->isCoId(id))
      leavePref[n] = coPref[thesolver->number(id)];

   left4X(n, id, 0, 1);
}

int SPxSteepPR::selectLeaveX(Real& best, int start, int incr)
{
   const Real* coPenalty_ptr = coPenalty.get_const_ptr();
   const Real* fTest = thesolver->fTest().get_const_ptr();
   //    const Real* low   = thesolver->lbBound();
   //    const Real* up    = thesolver->ubBound();
   const Real* p = leavePref.get_const_ptr();

   Real x;
   int selIdx;

   best = -infinity;
   selIdx = -1;

   //for(int i = 0; i < thesolver->dim() - start; i += incr)
 
   for (int i = thesolver->dim() - 1 - start; i >= 0; i -= incr)
   {
      x = fTest[i];

      if (x < -theeps)
      {         
         /**@todo this was an assert! is an assertion correct?*/
         // assert(coPenalty_ptr[i] >= theeps);
#ifndef NDEBUG
         if( coPenalty_ptr[i] < theeps )
         {
            std::cerr << "SPxSteepPR::selectLeaveX(): This shall not be!"
                      << std::endl;
            std::cerr << "i=" << i
                      << " x=" << x
                      << " coPenalty_ptr[i]=" << coPenalty_ptr[i]
                      << " theeps=" << theeps << std::endl;
            x = x * x / theeps * p[i];
         }
         else
#endif
            x = x * x / coPenalty_ptr[i] * p[i];
         if (x > best)
         {
            best = x;
            selIdx = i;
         }
      }
   }

   /*
       if(selIdx >= 0)
           std::cerr << fTest[selIdx] << std::endl;
    */ 
   return selIdx;
}

int SPxSteepPR::selectLeave()
{
   assert(isConsistent());
   Real best;

   lastIdx = selectLeaveX(best);
   if (lastIdx >= 0)
   {
      thesolver->basis().coSolve(thesolver->coPvec().delta(),
                                 thesolver->unitVector(lastIdx));
      workRhs.epsilon = accuracy;
      workRhs.setup_and_assign(thesolver->coPvec().delta());
      thesolver->setup4solve(&workVec, &workRhs);
   }

   /*
       if(thesolver->basis().iteration() < 10)
       {
           std::cerr.precision(16);
           std::cerr << lastIdx
                << '\t' << thesolver->fTest()[lastIdx]
                << '\t' << leavePref[lastIdx]
                << '\t' << coPenalty[lastIdx]
                << std::endl;
       }
    */

   return lastIdx;
}


/* Entering Simplex
 */
void SPxSteepPR::entered4X(
   SPxId, int n, int start2, int incr2, int start1, int incr1)
{
   assert(thesolver->type() == SoPlex::ENTER);

   if (n >= 0 && n < thesolver->dim())
   {
      Real delta = 2 + 1.0 / thesolver->basis().iteration();
      Real* coPenalty_ptr = coPenalty.get_ptr();
      Real* penalty_ptr = penalty.get_ptr();
      const Real* workVec_ptr = workVec.get_const_ptr();
      const Real* pVec = thesolver->pVec().delta().values();
      const IdxSet& pIdx = thesolver->pVec().idx();
      const Real* coPvec = thesolver->coPvec().delta().values();
      const IdxSet& coPidx = thesolver->coPvec().idx();
      Real xi_p = 1 / thesolver->fVec().delta()[n];
      int i, j;
      Real xi_ip, x;

      assert(thesolver->fVec().delta()[n] > thesolver->epsilon()
              || thesolver->fVec().delta()[n] < -thesolver->epsilon());

      for (j = coPidx.size() - 1 - start1; j >= 0; j -= incr1)
      {
         i = coPidx.index(j);
         xi_ip = xi_p * coPvec[i];
         x = coPenalty_ptr[i] += xi_ip * (xi_ip * pi_p - 2 * workVec_ptr[i]);
         /*
         if(x < 1)
             coPenalty_ptr[i] = 1 / (2-x);
         */
         if (x < delta)
            coPenalty_ptr[i] = delta;
         // coPenalty_ptr[i] = 1;
         else if (x > infinity)
            coPenalty_ptr[i] = 1 / thesolver->epsilon();
      }

      for (j = pIdx.size() - 1 - start2; j >= 0; j -= incr2)
      {
         i = pIdx.index(j);
         xi_ip = xi_p * pVec[i];
         x = penalty_ptr[i] += xi_ip * (xi_ip * pi_p
            - 2 * (thesolver->vector(i) * workVec));
         /*
         if(x < 1)
             penalty_ptr[i] = 1 / (2-x);
         */
         if (x < delta)
            penalty_ptr[i] = delta;
         // penalty_ptr[i] = 1;
         else if (x > infinity)
            penalty_ptr[i] = 1 / thesolver->epsilon();
      }
   }

   /*@
       if(thesolver->isId(id))
           penalty[   thesolver->number(id) ] *= 1.0001;
       else if(thesolver->isCoId(id))
           coPenalty[ thesolver->number(id) ] *= 1.0001;
   */
}

void SPxSteepPR::entered4(SPxId id, int n)
{
   entered4X(id, n, 0, 1, 0, 1);
}

SPxId SPxSteepPR::selectEnterX(
   Real& best, int start1, int incr1, int start2, int incr2)
{
   /*
       std::cerr << "selectEnter " << start1 << '(' << incr1 << ")\t"
                              << start2 << '(' << incr2 << ")\n";
    */
   const Real* p = pref.get_const_ptr();
   const Real* cp = coPref.get_const_ptr();
   const Real* test = thesolver->test().get_const_ptr();
   const Real* coTest = thesolver->coTest().get_const_ptr();
   const Real* penalty_ptr = penalty.get_const_ptr();
   const Real* coPenalty_ptr = coPenalty.get_const_ptr();
   Real x;
   int i, end;

   SPxId selId;
   best = -infinity;

   for (end = thesolver->coDim(), i = start2; i < end; i += incr2)
   {
      x = test[i];
      if (x < -theeps)
      {
         x *= x / penalty_ptr[i];
         x *= p[i];
         // x *= 1 + p[i];
         if (x > best)
         {
            best = x;
            selId = thesolver->id(i);
         }
      }
   }
   incr2 = end;

   for (end = thesolver->dim(), i = start1; i < end; i += incr1)
   {
      x = coTest[i];
      if (x < -theeps)
      {
         x *= x / coPenalty_ptr[i];
         x *= cp[i];
         // x *= 1 + cp[i];
         if (x > best)
         {
            best = x;
            selId = thesolver->coId(i);
         }
      }
   }

   assert(isConsistent());
   return selId;
}

SPxId SPxSteepPR::selectEnter()
{
   Real best;
   lastId = SPxSteepPR::selectEnterX(best);

   if (lastId.isValid())
   {
      SSVector& delta = thesolver->fVec().delta();

      thesolver->basis().solve4update(delta, thesolver->vector(lastId));

      // workRhs.epsilon = 0.1*accuracy;
      workRhs.epsilon = accuracy;
      workRhs.setup_and_assign(delta);
      pi_p = 1 + delta.length2();

      thesolver->setup4coSolve(&workVec, &workRhs);
   }

   return lastId;
}

void SPxSteepPR::addedVecs(int n)
{
   n = penalty.dim();
   pref.reSize (thesolver->coDim());
   penalty.reDim(thesolver->coDim());

   if (thesolver->type() == SoPlex::ENTER)
   {
      setupPrefs(thesolver->type());
      for (; n < penalty.dim(); ++n)
         penalty[n] = 2;
   }
   prefSetup = 0;
}

void SPxSteepPR::addedCoVecs(int n)
{
   n = coPenalty.dim();

   leavePref.reSize(thesolver->dim());
   coPref.reSize (thesolver->dim());
   setupPrefs(thesolver->type());

   workVec.reDim (thesolver->dim());
   coPenalty.reDim (thesolver->dim());
   for (; n < coPenalty.dim(); ++n)
      coPenalty[n] = 1;
   prefSetup = 0;
}

void SPxSteepPR::removedVec(int i)
{
   assert(thesolver != 0);
   penalty[i] = penalty[penalty.dim()];
   penalty.reDim(thesolver->coDim());
   prefSetup = 0;
}

void SPxSteepPR::removedVecs(const int perm[])
{
   assert(thesolver != 0);
   if (thesolver->type() == SoPlex::ENTER)
   {
      int i;
      int j = penalty.dim();
      for (i = 0; i < j; ++i)
         if (perm[i] >= 0)
            penalty[perm[i]] = penalty[i];
   }
   penalty.reDim(thesolver->coDim());
   prefSetup = 0;
}

void SPxSteepPR::removedCoVec(int i)
{
   assert(thesolver != 0);
   coPenalty[i] = coPenalty[coPenalty.dim()];
   coPenalty.reDim(thesolver->dim());
   prefSetup = 0;
}

void SPxSteepPR::removedCoVecs(const int perm[])
{
   assert(thesolver != 0);
   int i;
   int j = coPenalty.dim();
   for (i = 0; i < j; ++i)
      if (perm[i] >= 0)
         coPenalty[perm[i]] = coPenalty[i];
   coPenalty.reDim(thesolver->dim());
   prefSetup = 0;
}

bool SPxSteepPR::isConsistent() const
{
   if (thesolver != 0 && thesolver->type() == SoPlex::LEAVE && setup == EXACT)
   {
      int i;
      SSVector tmp(thesolver->dim(), thesolver->epsilon());
      Real x;
      for (i = thesolver->dim() - 1; i >= 0; --i)
      {
         thesolver->basis().coSolve(tmp, thesolver->unitVector(i));
         x = coPenalty[i] - tmp.length2();
         if (x > thesolver->delta() || -x > thesolver->delta())
         {
            std::cerr << "x[" << i << "] = " << x << std::endl;
         }
      }
   }

   if (thesolver != 0 && thesolver->type() == SoPlex::ENTER)
   {
      int i;
      for (i = thesolver->dim() - 1; i >= 0; --i)
         if (coPenalty[i] < thesolver->epsilon())
            return MSGinconsistent("SPxSteepPR");

      for (i = thesolver->coDim() - 1; i >= 0; --i)
         if (penalty[i] < thesolver->epsilon())
            return MSGinconsistent("SPxSteepPR");
   }
   return true;
}
} // namespace soplex

//-----------------------------------------------------------------------------
//Emacs Local Variables:
//Emacs mode:c++
//Emacs c-basic-offset:3
//Emacs tab-width:8
//Emacs indent-tabs-mode:nil
//Emacs End:
//-----------------------------------------------------------------------------

invivospxsteeppr__ZN6soplex10SPxSteepPR9entered4XENS_5SPxIdEiiiii_399

Real time: 3.997608e+11 -- Predicted time: 7.590169e+09

Callcount: 10000

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5539	### 1	12308.2959125675	6.166710e+05	6.474480e+05	4.7
5416	### 2	6088.38203954857	0.000000e+00	6.509444e+07	100.0

                /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1997-1999 Roland Wunderling                              */
/*                  1997-2002 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.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef SPEC_CPU
#pragma ident "@(#) $Id: spxsteeppr.cpp,v 1.18 2002/03/21 16:06:19 bzfkocht Exp $"
#endif

//#define DEBUGGING 1

#include <assert.h>
#include <iostream>

#include "spxdefines.h"
#include "message.h"
#include "spxsteeppr.h"
#include "random.h"

namespace soplex
{

#define EQ_PREF 1000

void SPxSteepPR::clear()
{
   thesolver = 0;
   prefSetup = 0;
}

void SPxSteepPR::load(SoPlex* base)
{
   thesolver = base;

   if (base)
   {
      workVec.clear();
      workVec.reDim(base->dim());
      workRhs.clear();
      workRhs.reDim(base->dim());

      leavePref.reSize(base->dim());
      coPref.reSize (base->dim());
      pref.reSize (base->coDim());
      prefSetup = 0;
   }
}

void SPxSteepPR::setType(SoPlex::Type type)
{
   int i;

   workRhs.epsilon = thesolver->epsilon();

   pref.reSize (thesolver->coDim());
   coPref.reSize(thesolver->dim());
   setupPrefs(type);

   if (setup == DEFAULT)
   {
      if (type == SoPlex::ENTER)
      {
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
            // coPenalty[i] = 10;
            coPenalty[i] = 2;
         penalty.reDim(thesolver->coDim());
         for (i = thesolver->coDim() - 1; i >= 0; --i)
            // penalty[i] = 10;
            penalty[i] = 1;
         // penalty[i] = 1 + thesolver->vector(i).size() / thesolver->dim();
      }
      else
      {
         assert(type == SoPlex::LEAVE);
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
         {
            // coPenalty[i] = 1;
            SPxId id = thesolver->basis().baseId(i);
            int n = thesolver->number(id);
            if (thesolver->isId(id))
               leavePref[i] = pref[n];
            else
               leavePref[i] = coPref[n];
            coPenalty[i] = 1 + thesolver->basis().baseVec(i).size()
                           / Real(thesolver->dim());
         }
      }
   }
   else
   {
      std::cerr << 
         "sorry, no exact setup for steepest edge multipliers implemented\n";

      if (type == SoPlex::ENTER)
      {
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
            coPenalty[i] = 1;
         penalty.reDim(thesolver->coDim());
         for (i = thesolver->coDim() - 1; i >= 0; --i)
            penalty[i] = 1 + thesolver->vector(i).length2();
      }
      else
      {
         assert(type == SoPlex::LEAVE);
         coPenalty.reDim(thesolver->dim());
         for (i = thesolver->dim() - 1; i >= 0; --i)
         {
            coPenalty[i] = 1 + thesolver->basis().baseVec(i).size()
                           / Real(thesolver->dim());
         }
      }
   }

   workVec.clear();
   workRhs.clear();
}

void SPxSteepPR::setupPrefs(Real mult, Real /*tie*/, Real /*cotie*/,
                             Real shift, Real coshift,
                             int rs, int cs, int re, int ce)
{
   Real *p;
   Real *cp;
   Real *end;
   // Real rtie, ctie;
   Real rshift, cshift;
   int i;

   if (thesolver->rep() == SoPlex::COLUMN)
   {
      cp = pref.get_ptr();
      p = coPref.get_ptr();
      // ctie = tie;
      // rtie = cotie;
      cshift = shift;
      rshift = coshift;
   }
   else
   {
      p = pref.get_ptr();
      cp = coPref.get_ptr();
      // rtie = tie;
      // ctie = cotie;
      rshift = shift;
      cshift = coshift;
   }

   if (re < 0)
      re = thesolver->nRows();
   for (i = re; --i >= rs;)
   {
      p[i] = rshift;
      //      p[i] += rtie * thesolver->rowVector(i).size() / Real(thesolver->nCols());
      //      p[i] += EQ_PREF * (thesolver->rhs(i) == thesolver->lhs(i));
      //      p[i] += EQ_PREF * (thesolver->rhs(i) >=  infinity
      //                     &&  thesolver->lhs(i) <= -infinity);
   }

   if (ce < 0)
      ce = thesolver->nCols();
   for (i = ce; --i >= cs;)
   {
      cp[i] = cshift;
      //      cp[i] += ctie * thesolver->colVector(i).size() / Real(thesolver->nRows());
      //      cp[i] += EQ_PREF * (thesolver->upper(i) == thesolver->lower(i));
      //      cp[i] += EQ_PREF * (thesolver->upper(i) >=  infinity
      //                      &&  thesolver->lower(i) <= -infinity);
   }


   i = 0;
   cp = coPref.get_ptr();
   end = cp + coPref.size();
   while (cp < end)
      *cp++ *= 1.0 - mult * i++;

   p = pref.get_ptr();
   end = p + pref.size();
   i = pref.size();
   while (p < end)
      *p++ *= 1.0 + mult * i--;
}

void SPxSteepPR::setupPrefs(SoPlex::Type tp)
{
   if (tp != prefSetup)
   {
      Real mult = 1e-8 / Real(1 + thesolver->dim() + thesolver->coDim());
      if (tp == SoPlex::ENTER)
         setupPrefs(-mult, -1e-5, -1e-5, 1, 1);
      else
         setupPrefs(mult, 1e-5, 1e-5, 1, 1);
      prefSetup = tp;
   }
}

void SPxSteepPR::setRep(SoPlex::Representation)
{
   if (workVec.dim() != thesolver->dim())
   {
      DVector tmp = penalty;
      penalty = coPenalty;
      coPenalty = tmp;

      workVec.clear();
      workVec.reDim(thesolver->dim());
   }
}

void SPxSteepPR::left4X(int n, SPxId id, int start, int incr)
{
   assert(thesolver->type() == SoPlex::LEAVE);

   if (id.isValid())
   {
      int i, j;
      Real x;
      // Real               delta         = 0.1;   // thesolver->epsilon();
      Real delta = 0.1 + 1.0 / thesolver->basis().iteration();
      Real* coPenalty_ptr = coPenalty.get_ptr();
      const Real* workVec_ptr = workVec.get_const_ptr();
      const Real* rhoVec = thesolver->fVec().delta().values();
      Real rhov_1 = 1 / rhoVec[n];
      Real beta_q = thesolver->coPvec().delta().length2()
                      * rhov_1 * rhov_1;

      assert(rhoVec[n] >= theeps || -rhoVec[n] >= theeps);

      //  Update #coPenalty# vector
      const IdxSet& rhoIdx = thesolver->fVec().idx();
      int len = thesolver->fVec().idx().size();
      for (i = len - 1 - start; i >= 0; i -= incr)
      {
         j = rhoIdx.index(i);
         x = coPenalty_ptr[j] += rhoVec[j]
            * (beta_q * rhoVec[j] - 2 * rhov_1 * workVec_ptr[j]);
         if (x < delta)
            // coPenalty_ptr[j] = delta / (1+delta-x);
            coPenalty_ptr[j] = delta;
         else if (x >= infinity)
            coPenalty_ptr[j] = 1 / theeps;
      }

      coPenalty_ptr[n] = beta_q;
      //@ coPenalty_ptr[n] = 0.999*beta_q;
      //@ coPenalty_ptr[n] = 1.001*beta_q;
   }

}

void SPxSteepPR::left4(int n, SPxId id)
{
   //  Update preference multiplier in #leavePref#
   if (thesolver->isId(id))
      leavePref[n] = pref[thesolver->number(id)];
   else if (thesolver->isCoId(id))
      leavePref[n] = coPref[thesolver->number(id)];

   left4X(n, id, 0, 1);
}

int SPxSteepPR::selectLeaveX(Real& best, int start, int incr)
{
   const Real* coPenalty_ptr = coPenalty.get_const_ptr();
   const Real* fTest = thesolver->fTest().get_const_ptr();
   //    const Real* low   = thesolver->lbBound();
   //    const Real* up    = thesolver->ubBound();
   const Real* p = leavePref.get_const_ptr();

   Real x;
   int selIdx;

   best = -infinity;
   selIdx = -1;

   //for(int i = 0; i < thesolver->dim() - start; i += incr)
 
   for (int i = thesolver->dim() - 1 - start; i >= 0; i -= incr)
   {
      x = fTest[i];

      if (x < -theeps)
      {         
         /**@todo this was an assert! is an assertion correct?*/
         // assert(coPenalty_ptr[i] >= theeps);
#ifndef NDEBUG
         if( coPenalty_ptr[i] < theeps )
         {
            std::cerr << "SPxSteepPR::selectLeaveX(): This shall not be!"
                      << std::endl;
            std::cerr << "i=" << i
                      << " x=" << x
                      << " coPenalty_ptr[i]=" << coPenalty_ptr[i]
                      << " theeps=" << theeps << std::endl;
            x = x * x / theeps * p[i];
         }
         else
#endif
            x = x * x / coPenalty_ptr[i] * p[i];
         if (x > best)
         {
            best = x;
            selIdx = i;
         }
      }
   }

   /*
       if(selIdx >= 0)
           std::cerr << fTest[selIdx] << std::endl;
    */ 
   return selIdx;
}

int SPxSteepPR::selectLeave()
{
   assert(isConsistent());
   Real best;

   lastIdx = selectLeaveX(best);
   if (lastIdx >= 0)
   {
      thesolver->basis().coSolve(thesolver->coPvec().delta(),
                                 thesolver->unitVector(lastIdx));
      workRhs.epsilon = accuracy;
      workRhs.setup_and_assign(thesolver->coPvec().delta());
      thesolver->setup4solve(&workVec, &workRhs);
   }

   /*
       if(thesolver->basis().iteration() < 10)
       {
           std::cerr.precision(16);
           std::cerr << lastIdx
                << '\t' << thesolver->fTest()[lastIdx]
                << '\t' << leavePref[lastIdx]
                << '\t' << coPenalty[lastIdx]
                << std::endl;
       }
    */

   return lastIdx;
}


/* Entering Simplex
 */
void SPxSteepPR::entered4X(
   SPxId, int n, int start2, int incr2, int start1, int incr1)
{
   assert(thesolver->type() == SoPlex::ENTER);

   if (n >= 0 && n < thesolver->dim())
   {
      Real delta = 2 + 1.0 / thesolver->basis().iteration();
      Real* coPenalty_ptr = coPenalty.get_ptr();
      Real* penalty_ptr = penalty.get_ptr();
      const Real* workVec_ptr = workVec.get_const_ptr();
      const Real* pVec = thesolver->pVec().delta().values();
      const IdxSet& pIdx = thesolver->pVec().idx();
      const Real* coPvec = thesolver->coPvec().delta().values();
      const IdxSet& coPidx = thesolver->coPvec().idx();
      Real xi_p = 1 / thesolver->fVec().delta()[n];
      int i, j;
      Real xi_ip, x;

      assert(thesolver->fVec().delta()[n] > thesolver->epsilon()
              || thesolver->fVec().delta()[n] < -thesolver->epsilon());

      for (j = coPidx.size() - 1 - start1; j >= 0; j -= incr1)
      {
         i = coPidx.index(j);
         xi_ip = xi_p * coPvec[i];
         x = coPenalty_ptr[i] += xi_ip * (xi_ip * pi_p - 2 * workVec_ptr[i]);
         /*
         if(x < 1)
             coPenalty_ptr[i] = 1 / (2-x);
         */
         if (x < delta)
            coPenalty_ptr[i] = delta;
         // coPenalty_ptr[i] = 1;
         else if (x > infinity)
            coPenalty_ptr[i] = 1 / thesolver->epsilon();
      }

      for (j = pIdx.size() - 1 - start2; j >= 0; j -= incr2)
      {
         i = pIdx.index(j);
         xi_ip = xi_p * pVec[i];
         x = penalty_ptr[i] += xi_ip * (xi_ip * pi_p
            - 2 * (thesolver->vector(i) * workVec));
         /*
         if(x < 1)
             penalty_ptr[i] = 1 / (2-x);
         */
         if (x < delta)
            penalty_ptr[i] = delta;
         // penalty_ptr[i] = 1;
         else if (x > infinity)
            penalty_ptr[i] = 1 / thesolver->epsilon();
      }
   }

   /*@
       if(thesolver->isId(id))
           penalty[   thesolver->number(id) ] *= 1.0001;
       else if(thesolver->isCoId(id))
           coPenalty[ thesolver->number(id) ] *= 1.0001;
   */
}

void SPxSteepPR::entered4(SPxId id, int n)
{
   entered4X(id, n, 0, 1, 0, 1);
}

SPxId SPxSteepPR::selectEnterX(
   Real& best, int start1, int incr1, int start2, int incr2)
{
   /*
       std::cerr << "selectEnter " << start1 << '(' << incr1 << ")\t"
                              << start2 << '(' << incr2 << ")\n";
    */
   const Real* p = pref.get_const_ptr();
   const Real* cp = coPref.get_const_ptr();
   const Real* test = thesolver->test().get_const_ptr();
   const Real* coTest = thesolver->coTest().get_const_ptr();
   const Real* penalty_ptr = penalty.get_const_ptr();
   const Real* coPenalty_ptr = coPenalty.get_const_ptr();
   Real x;
   int i, end;

   SPxId selId;
   best = -infinity;

   for (end = thesolver->coDim(), i = start2; i < end; i += incr2)
   {
      x = test[i];
      if (x < -theeps)
      {
         x *= x / penalty_ptr[i];
         x *= p[i];
         // x *= 1 + p[i];
         if (x > best)
         {
            best = x;
            selId = thesolver->id(i);
         }
      }
   }
   incr2 = end;

   for (end = thesolver->dim(), i = start1; i < end; i += incr1)
   {
      x = coTest[i];
      if (x < -theeps)
      {
         x *= x / coPenalty_ptr[i];
         x *= cp[i];
         // x *= 1 + cp[i];
         if (x > best)
         {
            best = x;
            selId = thesolver->coId(i);
         }
      }
   }

   assert(isConsistent());
   return selId;
}

SPxId SPxSteepPR::selectEnter()
{
   Real best;
   lastId = SPxSteepPR::selectEnterX(best);

   if (lastId.isValid())
   {
      SSVector& delta = thesolver->fVec().delta();

      thesolver->basis().solve4update(delta, thesolver->vector(lastId));

      // workRhs.epsilon = 0.1*accuracy;
      workRhs.epsilon = accuracy;
      workRhs.setup_and_assign(delta);
      pi_p = 1 + delta.length2();

      thesolver->setup4coSolve(&workVec, &workRhs);
   }

   return lastId;
}

void SPxSteepPR::addedVecs(int n)
{
   n = penalty.dim();
   pref.reSize (thesolver->coDim());
   penalty.reDim(thesolver->coDim());

   if (thesolver->type() == SoPlex::ENTER)
   {
      setupPrefs(thesolver->type());
      for (; n < penalty.dim(); ++n)
         penalty[n] = 2;
   }
   prefSetup = 0;
}

void SPxSteepPR::addedCoVecs(int n)
{
   n = coPenalty.dim();

   leavePref.reSize(thesolver->dim());
   coPref.reSize (thesolver->dim());
   setupPrefs(thesolver->type());

   workVec.reDim (thesolver->dim());
   coPenalty.reDim (thesolver->dim());
   for (; n < coPenalty.dim(); ++n)
      coPenalty[n] = 1;
   prefSetup = 0;
}

void SPxSteepPR::removedVec(int i)
{
   assert(thesolver != 0);
   penalty[i] = penalty[penalty.dim()];
   penalty.reDim(thesolver->coDim());
   prefSetup = 0;
}

void SPxSteepPR::removedVecs(const int perm[])
{
   assert(thesolver != 0);
   if (thesolver->type() == SoPlex::ENTER)
   {
      int i;
      int j = penalty.dim();
      for (i = 0; i < j; ++i)
         if (perm[i] >= 0)
            penalty[perm[i]] = penalty[i];
   }
   penalty.reDim(thesolver->coDim());
   prefSetup = 0;
}

void SPxSteepPR::removedCoVec(int i)
{
   assert(thesolver != 0);
   coPenalty[i] = coPenalty[coPenalty.dim()];
   coPenalty.reDim(thesolver->dim());
   prefSetup = 0;
}

void SPxSteepPR::removedCoVecs(const int perm[])
{
   assert(thesolver != 0);
   int i;
   int j = coPenalty.dim();
   for (i = 0; i < j; ++i)
      if (perm[i] >= 0)
         coPenalty[perm[i]] = coPenalty[i];
   coPenalty.reDim(thesolver->dim());
   prefSetup = 0;
}

bool SPxSteepPR::isConsistent() const
{
   if (thesolver != 0 && thesolver->type() == SoPlex::LEAVE && setup == EXACT)
   {
      int i;
      SSVector tmp(thesolver->dim(), thesolver->epsilon());
      Real x;
      for (i = thesolver->dim() - 1; i >= 0; --i)
      {
         thesolver->basis().coSolve(tmp, thesolver->unitVector(i));
         x = coPenalty[i] - tmp.length2();
         if (x > thesolver->delta() || -x > thesolver->delta())
         {
            std::cerr << "x[" << i << "] = " << x << std::endl;
         }
      }
   }

   if (thesolver != 0 && thesolver->type() == SoPlex::ENTER)
   {
      int i;
      for (i = thesolver->dim() - 1; i >= 0; --i)
         if (coPenalty[i] < thesolver->epsilon())
            return MSGinconsistent("SPxSteepPR");

      for (i = thesolver->coDim() - 1; i >= 0; --i)
         if (penalty[i] < thesolver->epsilon())
            return MSGinconsistent("SPxSteepPR");
   }
   return true;
}
} // namespace soplex

//-----------------------------------------------------------------------------
//Emacs Local Variables:
//Emacs mode:c++
//Emacs c-basic-offset:3
//Emacs tab-width:8
//Emacs indent-tabs-mode:nil
//Emacs End:
//-----------------------------------------------------------------------------

invivospxsteeppr__ZN6soplex10SPxSteepPR12selectEnterXERdiiii_449

Real time: 1.015031e+11 -- Predicted time: 1.070856e+11

Callcount: 10000

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
183	### 1	267.285634691307	6.495866e+06	6.439792e+06	0.8
311	### 2	145.703283345977	8.101373e+06	8.536292e+06	5.0
442	### 3	272.874139533207	8.850690e+06	9.543076e+06	7.2
548	### 4	550.54258921967	9.550274e+06	1.062053e+07	10.0
6953	### 5	1084.24621876226	1.048345e+07	1.114053e+07	5.8
3489	### 6	1613.17707702771	1.057515e+07	1.155054e+07	8.4
4214	### 7	1763.27600000269	1.130209e+07	1.189177e+07	4.9
5250	### 8	1695.62108551535	1.110851e+07	1.216278e+07	8.6
5257	### 9	1683.03717415	1.122964e+07	1.246985e+07	9.9
2346	### 10	907.678337658133	1.146369e+07	1.284838e+07	10.7

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivospxsteeppr_svectorZNK6soplex7SVectormlERKNS_6VectorE_293

Real time: 3.189254e+11 -- Predicted time: 1.094412e+11

Callcount: 696945308

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
137016610	### 1	81235485.3382647	1.830000e+02	7.200000e+01	60.6
499309737	### 2	67705517.5951541	1.120000e+02	8.000000e+01	28.5
252353915	### 3	152349787.673453	1.150000e+02	3.600000e+01	68.6
298264618	### 4	59288726.9620942	2.730000e+02	4.800000e+01	82.4
668498322	### 5	92428265.8257311	2.260000e+02	8.400000e+01	62.8
506529820	### 6	75543957.5131399	7.100000e+01	6.800000e+01	4.2
83408542	### 7	66536261.1235273	1.580000e+02	7.600000e+01	51.8
478473124	### 8	30755621.119756	2.540000e+02	2.880000e+02	11.8
227837634	### 9	73654457.8875532	1.150000e+02	3.200000e+01	72.1
148473471	### 10	576617.294972828	4.140000e+02	9.800000e+02	57.7

                /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1997-1999 Roland Wunderling                              */
/*                  1997-2002 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.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef SPEC_CPU
#pragma ident "@(#) $Id: ssvector.cpp,v 1.33 2005/01/12 12:00:03 bzfkocht Exp $"
#endif

#include <iostream>
#include <iomanip>
#include <assert.h>

#include "spxdefines.h"
#include "ssvector.h"
#include "svset.h"
#include "message.h"

/**@file ssvector.cpp
 * @todo There is a lot pointer arithmetic done here. It is not clear if
 *       this is an advantage at all. See all the function int() casts.
 * @todo Several operations like maxAbs could setup the vector while
 *       computing there result.
 */
namespace soplex
{
#define MARKER   1e-100

static const Real shortProductFactor = 0.5;

void SSVector::setMax(int newmax)
{
   assert(idx    != 0);
   assert(newmax != 0);
   assert(newmax >= IdxSet::size());

   // len = (newmax < IdxSet::max()) ? IdxSet::max() : newmax;
   len = newmax;

   spx_realloc(idx, len);
}

void SSVector::reDim (int newdim)
{
   for (int i = IdxSet::size() - 1; i >= 0; --i)
      if (index(i) >= newdim)
         remove(i);
   DVector::reDim(newdim);
   setMax(DVector::memSize() + 1);
   assert(isConsistent());
}

void SSVector::reMem(int newsize)
{
   DVector::reSize(newsize);
   assert(isConsistent());
   setMax(DVector::memSize() + 1);
}

void SSVector::clear ()
{
   if (isSetup())
   {
      for(int i = 0; i < num; ++i)
         val[idx[i]] = 0.0;
   }
   else
      Vector::clear();

   IdxSet::clear();
   setupStatus = true;
   assert(isConsistent());
}

void SSVector::setValue(int i, Real x)
{
   assert(i >= 0 && i < DVector::dim());

   if (isSetup())
   {
      int n = number(i);

      if (n < 0)
      {
         if (isNotZero(x, epsilon))
            IdxSet::add(1, &i);
      }
      else if (x == 0)
         clearNum(n);
   }
   val[i] = x;

   assert(isConsistent());
}

#ifdef USE_OLD // old version
void SSVector::setup()
{
   if (!isSetup())
   {
      IdxSet::clear();

      // #define      TWO_LOOPS
#ifdef  TWO_LOOPS
      int i = 0;
      int n = 0;
      int* id = idx;
      Real* v = val;
      const Real* end = val + dim();

      while (v < end)
      {
         id[n] = i++;
         n += (*v++ != 0);
      }

      Real x;
      int* ii = idx;
      int* last = idx + n;
      v = val;

      for (; id < last; ++id)
      {
         x = v[*id];
         if (isNotZero(x, epsilon))
            *ii++ = *id;
         else
            v[*id] = 0;
      }
      num = ii - idx;

#else

      if (dim() <= 1)
      {
         if (dim())
         {
            if (isNotZero(*val, epsilon))
               IdxSet::add(0);
            else
               *val = 0;
         }
      }
      else
      {
         int* ii = idx;
         Real* v = val;
         Real* end = v + dim() - 1;

         /* setze weissen Elefanten */
         Real last = *end;
         *end = MARKER;

         /* erstes element extra */
         if (isNotZero(*v, epsilon))
            *ii++ = 0;
         else
            *v = 0.0;

         for(;;)
         {
            while (*++v == 0.0)
               ;
            if (isNotZero(*v, epsilon))
            {
               *ii++ = int(v - val);
            }
            else
            {
               *v = 0.0;
               if (v == end)
                  break;
            }
         }
         /* fange weissen Elefanten wieder ein */
         if (isNotZero(last, epsilon))
         {
            *v = last;
            *ii++ = dim() - 1;
         }
         else
            *v = 0;

         num = int(ii - idx);
      }

#endif

      setupStatus = true;
      assert(isConsistent());
   }
}
#else // new version, not yet fully testet
void SSVector::setup()
{
   if (!isSetup())
   {
      IdxSet::clear();

      num = 0;

      for(int i = 0; i < dim(); ++i)
      {
         if (val[i] != 0.0)
         {
            if (isZero(val[i], epsilon))
               val[i] = 0.0;
            else
            {
               idx[num] = i;
               num++;
            }
         }
      }
      setupStatus = true;
      assert(isConsistent());
   }
}
#endif

SSVector& SSVector::operator+=(const Vector& vec)
{
   Vector::operator+=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

SSVector& SSVector::operator+=(const SVector& vec)
{
   Vector::operator+=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator+=(const SSVector& vec)
{
   for (int i = vec.size() - 1; i >= 0; --i)
      val[vec.index(i)] += vec.value(i);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#else
SSVector& SSVector::operator+=(const SSVector& vec)
{
   for (int i = 0; i < vec.size(); ++i)
      val[vec.index(i)] += vec.value(i);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#endif

SSVector& SSVector::operator-=(const Vector& vec)
{
   Vector::operator-=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

SSVector& SSVector::operator-=(const SVector& vec)
{
   Vector::operator-=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator-=(const SSVector& vec)
{
   if (vec.isSetup())
   {
      for (int i = vec.size() - 1; i >= 0; --i)
         val[vec.index(i)] -= vec.value(i);
   }
   else
   {
      Vector::operator-=(Vector(vec));
   }

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#else
SSVector& SSVector::operator-=(const SSVector& vec)
{
   if (vec.isSetup())
   {
      for (int i = 0; i < vec.size(); ++i)
         val[vec.index(i)] -= vec.value(i);
   }
   else
   {
      Vector::operator-=(Vector(vec));
   }

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#endif

#ifdef USE_OLD // old version
SSVector& SSVector::operator*=(Real x)
{
   assert(isSetup());

   for (int i = size() - 1; i >= 0; --i)
      val[index(i)] *= x;
   assert(isConsistent());
   return *this;
}
#else
SSVector& SSVector::operator*=(Real x)
{
   assert(isSetup());

   for (int i = 0; i < size(); ++i)
      val[index(i)] *= x;

   assert(isConsistent());

   return *this;
}
#endif

#ifdef USE_OLD // old
Real SSVector::maxAbs() const
{
   if (isSetup())
   {
      int* i = idx;
      int* end = idx + num;
      Real* v = val;
      Real absval = 0.0;

      for (; i < end; ++i)
      {
         Real x = v[*i];
         if (fabs(x) > absval)
            absval = fabs(x);
      }
      return absval;
   }
   else
      return Vector::maxAbs();
}
#else // new, not fully tested
Real SSVector::maxAbs() const
{
   if (isSetup())
   {
      Real maxabs = 0.0;

      for(int i = 0; i < num; ++i)
      {
         Real x = fabs(val[idx[i]]);

         if (x > maxabs)
            maxabs = x;
      }
      return maxabs;
   }
   else
      return Vector::maxAbs();
}
#endif // !0

Real SSVector::length2() const
{
   Real x = 0.0;

   if (isSetup())
   {
      for(int i = 0; i < num; ++i)
         x += val[idx[i]] * val[idx[i]];
   }
   else
      x = Vector::length2();

   return x;
}

Real SSVector::length() const
{
   return sqrt(length2());
}

#if 0 // buggy and not used
/* @todo check if really not used or if the Vector version is used instead.
 */
SSVector& SSVector::multAdd(Real xx, const SSVector& svec)
{
   if (svec.isSetup())
   {
      if (isSetup())
      {
         int i, j;
         Real x;

         for (i = svec.size() - 1; i >= 0; --i)
         {
            j = svec.index(i);
            if (val[j])
            {
               x = val[j] + xx * svec.value(i);
               if (isNotZero(x, epsilon))
                  val[j] = x;
               else
               {
                  val[j] = 0;
                  for (--i; i >= 0; --i)
                     val[svec.index(i)] += xx * svec.value(i);
                  unSetup();
                  break;
               }
            }
            else
            {
               x = xx * svec.value(i);
               if (isNotZero(x, epsilon))
               {
                  val[j] = x;
                  addIdx(j);
               }
            }
         }
      }
      else
         Vector::multAdd(xx, svec);
   }
   else
   {
      /**@todo this code does not work, because in is never something
       *       added to v. Also the idx will not be setup correctly
       *       Fortunately the whole function seems not to be called
       *       at all. 
       */
      assert(false);

      Real y;
      int* ii = idx;
      Real* v = val;
      Real* rv = static_cast<Real*>(svec.val);
      Real* last = rv + svec.dim() - 1;
      Real x = *last;

      *last = MARKER;
      for(;;)
      {
         while (!*rv)
         {
            ++rv;
            ++v;
         }
         y = *rv++ * xx;
         if (isNotZero(y, epsilon))
         {
            *ii++ = int(v - val);
            *v++ = y;
         }
         else if (rv == last)
            break;
         else
            v++;
      }
      *rv = x;

      x *= xx;
      if (isNotZero(x, epsilon))
      {
         *ii++ = int(v - val);
         *v = x;
      }
      num = int(ii - idx);

      setupStatus = true;
   }

   assert(isConsistent());
   return *this;
}
#endif // 0

/* @todo This function does not look good. MARKER is set but never really used.
 */
SSVector& SSVector::multAdd(Real xx, const SVector& svec)
{
   if (isSetup())
   {
      int i, j;
      Real x;
      Real* v = val;
      int adjust = 0;

      for (i = svec.size() - 1; i >= 0; --i)
      {
         j = svec.index(i);
         if (v[j])
         {
            x = v[j] + xx * svec.value(i);
            if (isNotZero(x, epsilon))
               v[j] = x;
            else
            {
               adjust = 1;
               v[j] = MARKER;
            }
         }
         else
         {
            x = xx * svec.value(i);
            if (isNotZero(x, epsilon))
            {
               v[j] = x;
               addIdx(j);
            }
         }
      }

      if (adjust)
      {
         int* iptr = idx;
         int* iiptr = idx;
         int* endptr = idx + num;
         for (; iptr < endptr; ++iptr)
         {
            x = v[*iptr];
            if (isNotZero(x, epsilon))
               *iiptr++ = *iptr;
            else
               v[*iptr] = 0;
         }
         num = int(iiptr - idx);
      }
   }
   else
      Vector::multAdd(xx, svec);

   assert(isConsistent());
   return *this;
}

SSVector& SSVector::multAdd(Real x, const Vector& vec)
{
   Vector::multAdd(x, vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator=(const SSVector& rhs)
{
   assert(rhs.isConsistent());

   if (this != &rhs)
   {
      clear();

      setMax(rhs.max());
      IdxSet::operator=(rhs);
      DVector::reDim(rhs.dim());

      if (rhs.isSetup())
      {
         for (int i = size() - 1; i >= 0; --i)
         {
            int j = index(i);
            val[j] = rhs.val[j];
         }
      }
      else
      {
         int* ii = idx;
         Real* v = val;
         Real* rv = static_cast<Real*>(rhs.val);
         Real* last = rv + rhs.dim() - 1;
         Real x = *last;
         
         *last = MARKER;
         for(;;)
         {
            while (!*rv)
            {
               ++rv;
               ++v;
            }
            if (isNotZero(*rv, epsilon))
            {
               *ii++ = int(v - val);
               *v++ = *rv++;
            }
            else if (rv == last)
               break;
            else
            {
               v++;
               rv++;
            }
         }
         *rv = x;
         
         if (isNotZero(x, epsilon))
         {
            *ii++ = int(v - val);
            *v++ = x;
         }
         num = int(ii - idx);
      }
      setupStatus = true;
   }
   assert(isConsistent());

   return *this;
}
#else // new version
SSVector& SSVector::operator=(const SSVector& rhs)
{
   assert(rhs.isConsistent());

   if (this != &rhs)
   {
      clear();
      epsilon = rhs.epsilon;
      setMax(rhs.max());
      DVector::reDim(rhs.dim());

      if (rhs.isSetup())
      {
         IdxSet::operator=(rhs);

         for(int i = 0; i < size(); ++i)
         {
            int j  = index(i);
            val[j] = rhs.val[j];
         }
      }
      else
      {
         num = 0;

         for(int i = 0; i < rhs.dim(); ++i)
         {
            if (isNotZero(rhs.val[i], epsilon))
            {
               val[i]       = rhs.val[i];
               idx[num]     = i;
               num++;
            }
         }
      }
      setupStatus = true;
   }
   assert(isConsistent());

   return *this;
}
#endif // 0

#ifdef USE_OLD // old version
void SSVector::setup_and_assign(SSVector& rhs)
{
   assert(rhs.isConsistent());

   clear();
   epsilon = rhs.epsilon;
   setMax(rhs.max());
   DVector::reDim(rhs.dim());

   if (rhs.isSetup())
   {
      int i, j;
      IdxSet::operator=(rhs);
      for (i = size() - 1; i >= 0; --i)
      {
         j = index(i);
         val[j] = rhs.val[j];
      }
   }
   else
   {
      int* ri = rhs.idx;
      int* ii = idx;
      Real* rv = rhs.val;
      Real* v = val;
      Real* last = rv + rhs.dim() - 1;
      Real x = *last;

      *last = MARKER;

      for(;;)
      {
         while (!*rv)
         {
            ++rv;
            ++v;
         }
         if (isNotZero(*rv, rhs.epsilon))
         {
            *ri++ = *ii++ = int(v - val);
            *v++ = *rv++;
         }
         else if (rv == last)
            break;
         else
         {
            v++;
            *rv++ = 0;
         }
      }

      if (isNotZero(x, rhs.epsilon))
      {
         *ri++ = *ii++ = int(v - val);
         *v++ = *rv = x;
      }
      else
         *rv = 0;
      num = rhs.num = int(ii - idx);
      rhs.setupStatus = true;
   }
   setupStatus = true;

   assert(isConsistent());
}
#else // new version
/* setup rhs and assign to this
 */
void SSVector::setup_and_assign(SSVector& rhs)
{
   clear();
   epsilon = rhs.epsilon;
   setMax(rhs.max());
   DVector::reDim(rhs.dim());

   if (rhs.isSetup())
   {
      IdxSet::operator=(rhs);
      
      for(int i = 0; i < size(); ++i)
      {
         int j  = index(i);
         val[j] = rhs.val[j];
      }
   }
   else
   {
      num = 0;

      for(int i = 0; i < rhs.dim(); ++i)
      {
         if (rhs.val[i] != 0.0)
         {
            if (isNotZero(rhs.val[i], epsilon))
            {
               rhs.idx[num] = i;
               idx[num]     = i;
               val[i]       = rhs.val[i];
               num++;
            }
            else
            {
               rhs.val[i] = 0.0;
            }
         }
      }
      rhs.num         = num;
      rhs.setupStatus = true;
   }
   setupStatus = true;

   assert(rhs.isConsistent());
   assert(isConsistent());
}
#endif // 0


SSVector& SSVector::operator=(const SVector& rhs)
{
   clear();
   return assign(rhs);
}

#ifdef USE_OLD // old version (buggy or optimization dependend)
SSVector& SSVector::assign(const SVector& rhs)
{
   assert(rhs.dim() <= Vector::dim());

   const SVector::Element* e = rhs.m_elem;
   int* p = idx;
   int i = rhs.size();

   while (i--)
   {
      val[*p = e->idx] = e->val;
      p += ((e++)->val != 0);
   }
   num = int(p - idx);
   setupStatus = true;

   assert(isConsistent());
   return *this;
}

#else // new version not yet fully testet, could be put into operator=()
SSVector& SSVector::assign(const SVector& rhs)
{
   assert(rhs.dim() <= Vector::dim());

   num = 0;

   for(int i = 0; i < rhs.size(); ++i)
   {
      int  k = rhs.index(i);
      Real v = rhs.value(i);

      if (isZero(v, epsilon))
         val[k] = 0.0;
      else
      {
         val[k]     = v;
         idx[num++] = k;
      }
   }
   setupStatus = true;

   assert(isConsistent());

   return *this;
}
#endif

SSVector& SSVector::assign2product1(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   const Real* vl = x.val;
   const int* xi = x.idx;

   int* ii = idx;
   SVector* svec = const_cast<SVector*>( & A[*xi] );
   const SVector::Element* e = &(svec->element(0));
   const SVector::Element* last = e + (num = svec->size());
   Real* v = val;
   Real y = vl[*xi];

   for (; e < last; ++e)
      v[ *ii++ = e->idx ] = y * e->val;

   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productShort(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   int i;
   int j;
   Real y;
   int* ii                      = idx;
   const Real* vl               = x.val;
   const int*  xi               = x.idx;
   SVector* svec                = const_cast<SVector*>( & A[*xi] );
   const SVector::Element* e    = &(svec->element(0));
   const SVector::Element* last = e + (num = svec->size());
   Real* v                      = val;
   Real xx                      = vl[*xi++];

   for (; e < last; ++e)
   {
      y      = xx * e->val;
      *ii    = e->idx;
      v[*ii] = y;
      ii    += (y != 0) ? 1 : 0;
   }

   for (i = x.size(); --i > 0;)
   {
      xx = vl[*xi];
      svec = const_cast<SVector*>( & A[*xi++] );
      e = &(svec->element(0));
      for (int k = svec->size(); --k >= 0;)
      {
         j   = e->idx;
         *ii = j;
         y   = v[j];
         ii  += (y == 0) ? 1 : 0;
         y   += xx * e->val;
         e++;
         v[j] = (y != 0) ? y : MARKER;
      }
   }

   int* is = idx;
   int* it = idx;
   for (; is < ii; ++is)
   {
      if (isNotZero(v[*is], epsilon))
         *it++ = *is;
      else
         v[*is] = 0;
   }
   num = int(it - idx);

   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productFull(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   int i;
   const Real* vl = x.val;
   const int* xi = x.idx;

   SVector* svec;
   const SVector::Element* elem;
   const SVector::Element* last;
   Real y;
   Real* v = val;

   for (i = x.size(); i-- > 0; ++xi)
   {
      svec = const_cast<SVector*>( & A[*xi] );
      elem = &(svec->element(0));
      last = elem + svec->size();
      y = vl[*xi];
      for (; elem < last; ++elem)
         v[elem->idx] += y * elem->val;
   }
   return *this;
}

SSVector& SSVector::assign2product4setup(const SVSet& A, const SSVector& x)
{
   assert(A.num() == x.dim());

   assert(x.isSetup());

   clear();

   if (x.size() == 1)
   {
      assign2product1(A, x);
      setupStatus = true;
   }

   else if (Real(x.size())*A.memSize() <= shortProductFactor*dim()*A.num()
             && isSetup())
   {
      assign2productShort(A, x);
      setupStatus = true;
   }

   else
   {
      assign2productFull(A, x);
      setupStatus = false;
   }
   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2product(const SSVector& x, const SVSet& A)
{
   assert(A.num() == dim());

   Real y;

   clear();

   for (int i = dim(); i-- > 0;)
   {
      y = A[i] * x;

      if (isNotZero(y, epsilon))
      {
         val[i] = y;
         IdxSet::addIdx(i);
      }
   }
   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productAndSetup(const SVSet& A, SSVector& x)
{
   if (x.isSetup())
      return assign2product4setup(A, x);

   SVector* svec;
   const SVector::Element* elem;
   const SVector::Element* last;
   Real y;
   Real* v = val;
   int* xi = x.idx;
   Real* xv = x.val;
   Real* end = xv + x.dim() - 1;

   /* setze weissen Elefanten */
   Real lastval = *end;
   *end = MARKER;

   for(;;)
   {
      while (!*xv)
         ++xv;
      if (isNotZero(*xv, epsilon))
      {
         y = *xv;
         svec = const_cast<SVector*>( & A[ *xi++ = int(xv - x.val) ] );
         elem = &(svec->element(0));
         last = elem + svec->size();
         for (; elem < last; ++elem)
            v[elem->idx] += y * elem->val;
      }
      else
      {
         *xv = 0;
         if (xv == end)
            break;
      }
      xv++;
   }

   /* fange weissen Elefanten wieder ein */
   if (isNotZero(lastval, epsilon))
   {
      y = *xv = lastval;
      svec = const_cast<SVector*>( & A[ *xi++ = int(xv - x.val) ] );
      elem = &(svec->element(0));
      last = elem + svec->size();
      for (; elem < last; ++elem)
         v[elem->idx] += y * elem->val;
   }
   else
      *xv = 0;

   x.num = int(xi - x.idx);
   x.setupStatus = true;
   setupStatus = false;

   assert(isConsistent());

   return *this;
}

bool SSVector::isConsistent() const
{
   if (Vector::dim() > IdxSet::max())
      return MSGinconsistent("SSVector");

   if (Vector::dim() < IdxSet::dim())
      return MSGinconsistent("SSVector");

   if (isSetup())
   {
      for (int i = 0; i < Vector::dim(); ++i)
      {
         int j = number(i);

         if (j < 0 && fabs(val[i]) > 0.0) 
         {
            std::cerr << "i= " << i << " idx= " << j << " val= " << std::setprecision(16) << val[i] << std::endl;
            return MSGinconsistent("SSVector");
         }
      }
   }
   return DVector::isConsistent() && IdxSet::isConsistent();
}
} // namespace soplex

//-----------------------------------------------------------------------------
//Emacs Local Variables:
//Emacs mode:c++
//Emacs c-basic-offset:3
//Emacs tab-width:8
//Emacs indent-tabs-mode:nil
//Emacs End:
//-----------------------------------------------------------------------------

invivossvector__ZN6soplex8SSVector5clearEv_74

Real time: 1.639992e+10 -- Predicted time: 1.690773e+10

Callcount: 90023

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
64495	### 1	14398709.8408546	1.160000e+02	2.400000e+01	79.3
47902	### 2	6079.15679087496	2.506512e+06	2.637984e+06	4.9

                /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1997-1999 Roland Wunderling                              */
/*                  1997-2002 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.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef SPEC_CPU
#pragma ident "@(#) $Id: enter.cpp,v 1.16 2002/03/21 16:06:17 bzfkocht Exp $"
#endif

// #define DEBUGGING 1

/*      \SubSection{Updating the Basis for Entering Variables}
 */
#include <assert.h>

#include "spxdefines.h"
#include "soplex.h"
#include "spxratiotester.h"

namespace soplex
{

/*
In the entering simplex algorithms (i.e. iteratively a vector is selected to
{\em enter} the simplex basis as in the dual rowwise and primal columnwise case)
let $A$ denote the current basis, $x$ and entering vector and $f$ the
feasibility vector. For a feasible basis $l \le f \le u$ holds.  For the rowwise
case $f$ is obtained by solving $f^T = c^T A^{-1}$, wherease in columnwisecase
$f = A^{-1} b$.
 
Let us further consider the rowwise case. Exchanging $x$ with the $i$-th
vector of $A$ yields
\begin{equation}\label{update.eq}
    A^{(i)} = E_i A \hbox{, with } E_i = I + e_i (x^T A^{-1} - e_i^T).
\end{equation}
With $E_i^{-1} = I + e_i \frac{e_i^T - \delta^T}{\delta}$, $\delta^T = x^T A^{-1}$ one
gets the new feasibility vector
\begin{eqnarray*}
        (f^{(i)})^T
    &=& c^T (A^{(i)})^{-1}      \\
    &=& c^T A^{-1} + c^T A^{-1} e_i \frac{e_i^T - \delta^T}{\delta_i} \\
    &=& f^T + \frac{f_i}{\delta_i} e_i^T - \frac{f_i}{\delta_i} \delta^T.       \\
\end{eqnarray*}
The selection of the leaving vector $i^*$ for the basis must ensure, that for
all $j \ne i^*$ $f^{(i^*)}_j$ remains within its bounds $l_j$ and $u_j$.
 */


/*
    Testing all values of |pVec| against its bounds. If $i$, say, is violated
    the violation is saved as negative value in |theTest[i]|.
 */
Real SoPlex::test(int i, SPxBasis::Desc::Status stat) const
{
   METHOD( "SoPlex::test()" );
   assert(type() == ENTER);
   assert(!isBasic(stat));

   Real x;

   switch (stat)
   {
   case SPxBasis::Desc::D_FREE:
   case SPxBasis::Desc::D_ON_BOTH:
      assert(rep() == ROW);
      x = (*thePvec)[i] - lhs(i);
      if (x < 0)
         return x;
      // no break: next is else case
   case SPxBasis::Desc::D_ON_LOWER:
      assert(rep() == ROW);
      return rhs(i) - (*thePvec)[i];
   case SPxBasis::Desc::D_ON_UPPER:
      assert(rep() == ROW);
      return (*thePvec)[i] - lhs(i);

   case SPxBasis::Desc::P_ON_UPPER:
      assert(rep() == COLUMN);
      return maxObj(i) - (*thePvec)[i];
   case SPxBasis::Desc::P_ON_LOWER:
      assert(rep() == COLUMN);
      return (*thePvec)[i] - maxObj(i);
   case SPxBasis::Desc::P_FREE :
      x = maxObj(i) - (*thePvec)[i];
      return (x < 0) ? x : -x;

   default:
      return 0;
   }
}

void SoPlex::computeTest()
{
   METHOD( "SoPlex::computeTest()" );
   int i;
   const SPxBasis::Desc& ds = desc();

   for (i = coDim() - 1; i >= 0; --i)
   {
      SPxBasis::Desc::Status stat = ds.status(i);
      if (isBasic(stat))
         theTest[i] = 0;
      else
         theTest[i] = test(i, stat);
   }
}

Real SoPlex::computePvec(int i)
{
   METHOD( "SoPlex::computePvec()" );
   return (*thePvec)[i] = vector(i) * (*theCoPvec);
}

Real SoPlex::computeTest(int i)
{
   METHOD( "SoPlex::computeTest()" );
   SPxBasis::Desc::Status stat = desc().status(i);
   if (isBasic(stat))
      return theTest[i] = 0;
   else
      return theTest[i] = test(i, stat);
}

/*
    Testing all values of #coPvec# against its bounds. If $i$, say, is violated
    the violation is saved as negative value in |theCoTest[i]|.
 */
Real SoPlex::coTest(int i, SPxBasis::Desc::Status stat) const
{
   METHOD( "SoPlex::coTest()" );
   assert(type() == ENTER);
   assert(!isBasic(stat));

   Real x;

   switch (stat)
   {
   case SPxBasis::Desc::D_FREE:
   case SPxBasis::Desc::D_ON_BOTH :
      assert(rep() == ROW);
      x = (*theCoPvec)[i] - SPxLP::lower(i);
      if (x < 0)
         return x;
      // no break: next is else case
   case SPxBasis::Desc::D_ON_LOWER:
      assert(rep() == ROW);
      return SPxLP::upper(i) - (*theCoPvec)[i];
   case SPxBasis::Desc::D_ON_UPPER:
      assert(rep() == ROW);
      return (*theCoPvec)[i] - SPxLP::lower(i);

   case SPxBasis::Desc::P_ON_UPPER:
      assert(rep() == COLUMN);
      return (*theCoPvec)[i] - 0;             // slacks !
   case SPxBasis::Desc::P_ON_LOWER:
      assert(rep() == COLUMN);
      return 0 - (*theCoPvec)[i];             // slacks !

   default:
      return 0;
   }
}

void SoPlex::computeCoTest()
{
   METHOD( "SoPlex::computeCoTest()" );
   int i;
   const SPxBasis::Desc& ds = desc();

   for (i = dim() - 1; i >= 0; --i)
   {
      SPxBasis::Desc::Status stat = ds.coStatus(i);
      if (isBasic(stat))
         theCoTest[i] = 0;
      else
         theCoTest[i] = coTest(i, stat);
   }
}


/*
    The following methods require propersy initialized vectors |fVec| and
    #coPvec#.
 */
void SoPlex::updateTest()
{
   METHOD( "SoPlex::updateTest()" );
   thePvec->delta().setup();

   const IdxSet& idx = thePvec->idx();
   const SPxBasis::Desc& ds = desc();

   int i;
   for (i = idx.size() - 1; i >= 0; --i)
   {
      int j = idx.index(i);
      SPxBasis::Desc::Status stat = ds.status(j);
      if (!isBasic(stat))
         theTest[j] = test(j, stat);
      else
         theTest[j] = 0;
   }
}

void SoPlex::updateCoTest()
{
   METHOD( "SoPlex::updateCoTest()" );
   theCoPvec->delta().setup();

   const IdxSet& idx = theCoPvec->idx();
   const SPxBasis::Desc& ds = desc();

   int i;
   for (i = idx.size() - 1; i >= 0; --i)
   {
      int j = idx.index(i);
      SPxBasis::Desc::Status stat = ds.coStatus(j);
      if (!isBasic(stat))
         theCoTest[j] = coTest(j, stat);
      else
         theCoTest[j] = 0;
   }
}



/*  \Section{Compute statistics on entering variable}
    Here is a list of variables relevant when including |Id| to the basis.
    They are computed by |computeEnterStats()|.
 */
void SoPlex::getEnterVals
(
   SPxId enterId,
   Real& enterTest,
   Real& enterUB,
   Real& enterLB,
   Real& enterVal,
   Real& enterMax,
   Real& enterPric,
   SPxBasis::Desc::Status& enterStat,
   Real& enterRO
)
{
   METHOD( "SoPlex::getEnterVals()" );
   int enterIdx;
   SPxBasis::Desc& ds = desc();

   if (enterId.isSPxColId())
   {
      enterIdx = number(SPxColId(enterId));
      enterStat = ds.colStatus(enterIdx);
      assert(!isBasic(enterStat));

      /*      For an #Id# to enter the basis we better recompute the Test value.
       */
      if (rep() == COLUMN)
      {
         computePvec(enterIdx);
         enterTest = computeTest(enterIdx);
         theTest[enterIdx] = 0;
      }
      else
      {
         enterTest = coTest()[enterIdx];
         theCoTest[enterIdx] = 0;
      }

      switch (enterStat)
      {
         // primal/columnwise cases:
      case SPxBasis::Desc::P_ON_UPPER :
         assert( rep() == COLUMN );
         enterUB = theUCbound[enterIdx];
         enterLB = theLCbound[enterIdx];
         enterVal = enterUB;
         enterMax = enterLB - enterUB;
         enterPric = (*thePvec)[enterIdx];
         enterRO = maxObj(enterIdx);
         if( enterLB <= -infinity )
            ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_LOWER;
         else if( EQ( enterLB, enterUB ) )
            ds.colStatus(enterIdx) = SPxBasis::Desc::D_FREE;
         else
            ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
         break;
      case SPxBasis::Desc::P_ON_LOWER :
         assert( rep() == COLUMN );
         enterUB = theUCbound[enterIdx];
         enterLB = theLCbound[enterIdx];
         enterVal = enterLB;
         enterMax = enterUB - enterLB;
         enterPric = (*thePvec)[enterIdx];
         enterRO = maxObj(enterIdx);
         if( enterUB >= infinity )
            ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_UPPER;
         else if( EQ( enterLB, enterUB ) )
            ds.colStatus(enterIdx) = SPxBasis::Desc::D_FREE;
         else
            ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
         break;
      case SPxBasis::Desc::P_FREE :
         assert( rep() == COLUMN );
         enterUB = theUCbound[enterIdx];
         enterLB = theLCbound[enterIdx];
         enterVal = 0;
         enterPric = (*thePvec)[enterIdx];
         enterRO = maxObj(enterIdx);
         ds.colStatus(enterIdx) = SPxBasis::Desc::D_UNDEFINED;
         enterMax = (enterRO > 0) ? infinity : -infinity;
         break;

         // dual/rowwise cases:
      case SPxBasis::Desc::D_ON_UPPER :
         assert( rep() == ROW );
         assert(theUCbound[enterIdx] < infinity);
         enterUB = theUCbound[enterIdx];
         enterLB = -infinity;
         enterMax = -infinity;
         enterVal = enterUB;
         enterPric = (*theCoPvec)[enterIdx];
         enterRO = SPxLP::lower(enterIdx);
         ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
         break;
      case SPxBasis::Desc::D_ON_LOWER :
         assert( rep() == ROW );
         assert(theLCbound[enterIdx] > -infinity);
         enterLB = theLCbound[enterIdx];
         enterUB = infinity;
         enterMax = infinity;
         enterVal = enterLB;
         enterPric = (*theCoPvec)[enterIdx];
         enterRO = SPxLP::upper(enterIdx);
         ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
         break;
      case SPxBasis::Desc::D_FREE:
         assert( rep() == ROW );
         assert(SPxLP::lower(enterIdx) == SPxLP::upper(enterIdx));
         enterUB = infinity;
         enterLB = -infinity;
         enterVal = 0;
         enterRO = SPxLP::upper(enterIdx);
         enterPric = (*theCoPvec)[enterIdx];
         if (enterPric > enterRO)
            enterMax = infinity;
         else
            enterMax = -infinity;
         ds.colStatus(enterIdx) = SPxBasis::Desc::P_FIXED;
         break;
      case SPxBasis::Desc::D_ON_BOTH :
         assert( rep() == ROW );
         enterPric = (*theCoPvec)[enterIdx];
         if (enterPric > SPxLP::upper(enterIdx))
         {
            enterLB = theLCbound[enterIdx];
            enterUB = infinity;
            enterMax = infinity;
            enterVal = enterLB;
            enterRO = SPxLP::upper(enterIdx);
            ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
         }
         else
         {
            enterUB = theUCbound[enterIdx];
            enterVal = enterUB;
            enterRO = SPxLP::lower(enterIdx);
            enterLB = -infinity;
            enterMax = -infinity;
            ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
         }
         break;
      default:
         abort();
      }
      DEBUG({ std::cout << "SoPlex::getEnterVals() : col " << enterIdx
                        << ": " << enterStat
                        << " -> " << ds.colStatus(enterIdx)
                        << std::endl; });
   }

   else
   {
      assert(enterId.isSPxRowId());
      enterIdx = number(SPxRowId(enterId));
      enterStat = ds.rowStatus(enterIdx);
      assert(!isBasic(enterStat));

      /*      For an #Id# to enter the basis we better recompute the Test value.
       */
      if (rep() == ROW)
      {
         computePvec(enterIdx);
         enterTest = computeTest(enterIdx);
         theTest[enterIdx] = 0;
      }
      else
      {
         enterTest = coTest()[enterIdx];
         theCoTest[enterIdx] = 0;
      }

      switch (enterStat)
      {
         // primal/columnwise cases:
      case SPxBasis::Desc::P_ON_UPPER :
         assert( rep() == COLUMN );
         enterUB = theURbound[enterIdx];
         enterLB = theLRbound[enterIdx];
         enterVal = enterLB;
         enterMax = enterUB - enterLB;
         enterPric = (*theCoPvec)[enterIdx];
         enterRO = 0;
         if( enterUB >= infinity )
            ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_LOWER;
         else if( EQ( enterLB, enterUB ) )
            ds.rowStatus(enterIdx) = SPxBasis::Desc::D_FREE;
         else
            ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
         break;
      case SPxBasis::Desc::P_ON_LOWER :
         assert( rep() == COLUMN );
         enterUB = theURbound[enterIdx];
         enterLB = theLRbound[enterIdx];
         enterVal = enterUB;
         enterMax = enterLB - enterUB;
         enterPric = (*theCoPvec)[enterIdx];
         enterRO = 0;
         if( enterLB <= -infinity )
            ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_UPPER;
         else if( EQ( enterLB, enterUB ) )
            ds.rowStatus(enterIdx) = SPxBasis::Desc::D_FREE;
         else
            ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
         break;
      case SPxBasis::Desc::P_FREE :
         assert( rep() == COLUMN );
#if 1
         abort();
#else
         std::cerr << __FILE__ << __LINE__ << "ERROR: not yet debugged!\n";
         enterPric = (*theCoPvec)[enterIdx];
         enterRO = 0;
         ds.rowStatus(enterIdx) = SPxBasis::Desc::D_UNDEFINED;
#endif
         break;

         // dual/rowwise cases:
      case SPxBasis::Desc::D_ON_UPPER :
         assert( rep() == ROW );
         assert(theURbound[enterIdx] < infinity);
         enterUB = theURbound[enterIdx];
         enterLB = -infinity;
         enterVal = enterUB;
         enterMax = -infinity;
         enterPric = (*thePvec)[enterIdx];
         enterRO = lhs(enterIdx);
         ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
         break;
      case SPxBasis::Desc::D_ON_LOWER :
         assert( rep() == ROW );
         assert(theLRbound[enterIdx] > -infinity);
         enterLB = theLRbound[enterIdx];
         enterUB = infinity;
         enterVal = enterLB;
         enterMax = infinity;
         enterPric = (*thePvec)[enterIdx];
         enterRO = rhs(enterIdx);
         ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
         break;
      case SPxBasis::Desc::D_FREE:
         assert( rep() == ROW );
         assert(rhs(enterIdx) == lhs(enterIdx));
         enterUB = infinity;
         enterLB = -infinity;
         enterVal = 0;
         enterPric = (*thePvec)[enterIdx];
         enterRO = rhs(enterIdx);
         enterMax = (enterPric > enterRO) ? infinity : -infinity;
         ds.rowStatus(enterIdx) = SPxBasis::Desc::P_FIXED;
         break;
      case SPxBasis::Desc::D_ON_BOTH :
         assert( rep() == ROW );
         enterPric = (*thePvec)[enterIdx];
         if (enterPric > rhs(enterIdx))
         {
            enterLB = theLRbound[enterIdx];
            enterVal = enterLB;
            enterUB = infinity;
            enterMax = infinity;
            enterRO = rhs(enterIdx);
            ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
         }
         else
         {
            enterUB = theURbound[enterIdx];
            enterVal = enterUB;
            enterLB = -infinity;
            enterMax = -infinity;
            enterRO = lhs(enterIdx);
            ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
         }
         break;

      default:
         abort();
      }
      DEBUG({ std::cout << "SoPlex::getEnterVals() : row " << enterIdx
                        << ": " << enterStat
                        << " -> " << ds.rowStatus(enterIdx)
                        << std::endl; });
   }
}

/*      process leaving variable
 */
void SoPlex::getEnterVals2
(
   int leaveIdx,
   Real enterMax,
   Real& leavebound
)
{
   METHOD( "SoPlex::getEnterVals2()" );
   int idx;
   SPxBasis::Desc& ds = desc();
   SPxId leftId = baseId(leaveIdx);

   if (leftId.isSPxRowId())
   {
      idx = number(SPxRowId(leftId));
      SPxBasis::Desc::Status leaveStat = ds.rowStatus(idx);

      switch (leaveStat)
      {
      case SPxBasis::Desc::P_FIXED :
         assert(rep() == ROW);
         abort();
      case SPxBasis::Desc::P_ON_UPPER :
         assert(rep() == ROW);
         leavebound = theLBbound[leaveIdx];
         theLRbound[idx] = leavebound;
         ds.rowStatus(idx) = dualRowStatus(idx);
         break;
      case SPxBasis::Desc::P_ON_LOWER :
         assert(rep() == ROW);
         leavebound = theUBbound[leaveIdx];
         theURbound[idx] = leavebound;
         ds.rowStatus(idx) = dualRowStatus(idx);
         break;
      case SPxBasis::Desc::P_FREE :
         assert(rep() == ROW);
#if 1
         abort();
#else
         std::cerr << __FILE__ << __LINE__ << "ERROR: not yet debugged!\n";
         if ((*theCoPvec)[leaveIdx] - theLBbound[leaveIdx] <
              theUBbound[leaveIdx] - (*theCoPvec)[leaveIdx])
         {
            leavebound = theLBbound[leaveIdx];
            theLRbound[idx] = leavebound;
         }
         else
         {
            leavebound = theUBbound[leaveIdx];
            theURbound[idx] = leavebound;
         }
         ds.rowStatus(idx) = SPxBasis::Desc::D_UNDEFINED;
         break;
#endif

         // primal/columnwise cases:
      case SPxBasis::Desc::D_UNDEFINED :
         assert(rep() == COLUMN);
         abort();
      case SPxBasis::Desc::D_FREE :
         assert(rep() == COLUMN);
         if (theFvec->delta()[leaveIdx] * enterMax < 0)
            leavebound = theUBbound[leaveIdx];
         else
            leavebound = theLBbound[leaveIdx];
         theLRbound[idx] = leavebound;
         theURbound[idx] = leavebound;
         ds.rowStatus(idx) = SPxBasis::Desc::P_FIXED;
         break;
      case SPxBasis::Desc::D_ON_UPPER :
         assert(rep() == COLUMN);
         leavebound = theUBbound[leaveIdx];
         theURbound[idx] = leavebound;
         ds.rowStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
         break;
      case SPxBasis::Desc::D_ON_LOWER :
         assert(rep() == COLUMN);
         leavebound = theLBbound[leaveIdx];
         theLRbound[idx] = leavebound;
         ds.rowStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
         break;
      case SPxBasis::Desc::D_ON_BOTH :
         assert(rep() == COLUMN);
         if (enterMax * theFvec->delta()[leaveIdx] < 0)
         {
            leavebound = theUBbound[leaveIdx];
            theURbound[idx] = leavebound;
            ds.rowStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
         }
         else
         {
            leavebound = theLBbound[leaveIdx];
            theLRbound[idx] = leavebound;
            ds.rowStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
         }
         break;

      default:
         abort();
      }
      DEBUG({ std::cout << "SoPlex::getEnterVals2(): row " << idx
                        << ": " << leaveStat
                        << " -> " << ds.rowStatus(idx)
                        << std::endl; });
   }

   else
   {
      assert(leftId.isSPxColId());
      idx = number(SPxColId(leftId));
      SPxBasis::Desc::Status leaveStat = ds.colStatus(idx);

      switch (leaveStat)
      {
      case SPxBasis::Desc::P_ON_UPPER :
         assert(rep() == ROW);
         leavebound = theLBbound[leaveIdx];
         theLCbound[idx] = leavebound;
         ds.colStatus(idx) = dualColStatus(idx);
         break;
      case SPxBasis::Desc::P_ON_LOWER :
         assert(rep() == ROW);
         leavebound = theUBbound[leaveIdx];
         theUCbound[idx] = leavebound;
         ds.colStatus(idx) = dualColStatus(idx);
         break;
      case SPxBasis::Desc::P_FREE :
         assert(rep() == ROW);
         if (theFvec->delta()[leaveIdx] * enterMax > 0)
         {
            leavebound = theLBbound[leaveIdx];
            theLCbound[idx] = leavebound;
         }
         else
         {
            leavebound = theUBbound[leaveIdx];
            theUCbound[idx] = leavebound;
         }
         ds.colStatus(idx) = SPxBasis::Desc::D_UNDEFINED;
         break;
      case SPxBasis::Desc::P_FIXED:
         assert(rep() == ROW);
         abort();

         // primal/columnwise cases:
      case SPxBasis::Desc::D_FREE :
         assert(rep() == COLUMN);
         if (theFvec->delta()[leaveIdx] * enterMax > 0)
            leavebound = theLBbound[leaveIdx];
         else
            leavebound = theUBbound[leaveIdx];
         theUCbound[idx] =
            theLCbound[idx] = leavebound;
         ds.colStatus(idx) = SPxBasis::Desc::P_FIXED;
         break;
      case SPxBasis::Desc::D_ON_UPPER :
         assert(rep() == COLUMN);
         leavebound = theLBbound[leaveIdx];
         theLCbound[idx] = leavebound;
         ds.colStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
         break;
      case SPxBasis::Desc::D_ON_LOWER :
         assert(rep() == COLUMN);
         leavebound = theUBbound[leaveIdx];
         theUCbound[idx] = leavebound;
         ds.colStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
         break;
      case SPxBasis::Desc::D_ON_BOTH :
      case SPxBasis::Desc::D_UNDEFINED :
         assert(rep() == COLUMN);
         if (enterMax * theFvec->delta()[leaveIdx] < 0)
         {
            leavebound = theUBbound[leaveIdx];
            theUCbound[idx] = leavebound;
            ds.colStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
         }
         else
         {
            leavebound = theLBbound[leaveIdx];
            theLCbound[idx] = leavebound;
            ds.colStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
         }
         break;

      default:
         abort();
      }
      DEBUG({ std::cout << "SoPlex::getEnterVals2(): col " << idx
                        << ": " << leaveStat
                        << " -> " << ds.colStatus(idx)
                        << std::endl; });
   }
}


void
SoPlex::ungetEnterVal(
   SPxId enterId,
   SPxBasis::Desc::Status enterStat,
   Real leaveVal,
   const SVector& vec
)
{
   METHOD( "SoPlex::ungetEnterVal()" );
   int enterIdx;
   SPxBasis::Desc& ds = desc();

   if (enterId.isSPxColId())
   {
      enterIdx = number(SPxColId(enterId));
      if (enterStat == SPxBasis::Desc::P_ON_UPPER)
         ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
      else
         ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
      theFrhs->multAdd(leaveVal, vec);
   }
   else
   {
      enterIdx = number(SPxRowId(enterId));
      assert(enterId.isSPxRowId());
      if (enterStat == SPxBasis::Desc::P_ON_UPPER)
         ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
      else
         ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
      (*theFrhs)[enterIdx] += leaveVal;
   }
   if (isId(enterId))
      theTest[enterIdx] = 0;
   else
      theCoTest[enterIdx] = 0;
}

void SoPlex::rejectEnter(
   SPxId enterId,
   Real enterTest,
   SPxBasis::Desc::Status enterStat
)
{
   METHOD( "SoPlex::rejectEnter()" );
   int enterIdx = number(enterId);
   if (isId(enterId))
   {
      theTest[enterIdx] = enterTest;
      desc().status(enterIdx) = enterStat;
   }
   else
   {
      theCoTest[enterIdx] = enterTest;
      desc().coStatus(enterIdx) = enterStat;
   }
}

int SoPlex::enter(SPxId& enterId)
{
   METHOD( "SoPlex::enter()" );
   assert(enterId.isValid());
   assert(type() == ENTER);
   assert(initialized);

   Real enterTest;      // correct test value of entering var
   Real enterUB;        // upper bound of entering variable
   Real enterLB;        // lower bound of entering variable
   Real enterVal;       // current value of entering variable
   Real enterMax;       // maximum value for entering shift
   Real enterPric;      // priced value of entering variable
   SPxBasis::Desc::Status enterStat;      // status of entering variable
   Real enterRO;        // rhs/obj of entering variable
   const SVector* enterVec = enterVector(enterId);

   getEnterVals(enterId, enterTest, enterUB, enterLB,
                enterVal, enterMax, enterPric, enterStat, enterRO);
   //@ if(enterTest > ((theShift>epsilon()) ? -delta() : -epsilon()))
   if (enterTest > -epsilon())
   {
      rejectEnter(enterId, enterTest, enterStat);
      change(-1, enterId, enterVec);
      DEBUG( std::cerr << "rejecting false enter pivot" << std::endl; );
      return 1;
   }

   /*  Before performing the actual basis update, we must determine, how this
       is to be accomplished.
    */
   if (theFvec->delta().isSetup() && theFvec->delta().size() == 0)
      SPxBasis::solve4update(theFvec->delta(), *enterVec);
#ifndef NDEBUG
   else
   {
      DVector tmp(dim());
      tmp = reinterpret_cast<DVector&>(theFvec->delta());
      multBaseWith(tmp);
      tmp -= *enterVec;
      if (tmp.length() > delta())
         std::cerr << "fVec updated error = " << tmp.length() << std::endl;
   }
#endif  // NDEBUG

   if (m_numCycle > m_maxCycle)
   {
      if (-enterMax > 0)
         perturbMaxEnter();
      else
         perturbMinEnter();
   }

   Real leaveVal = -enterMax;
   int leaveIdx = theratiotester->selectLeave(leaveVal);


   /*
       We now tried to find a variable to leave the basis. If one has been
       found, a regular basis update is to be performed.
    */
   if (leaveIdx >= 0)
   {
      if (leaveVal < delta() && leaveVal > -delta())
         m_numCycle += (theUBbound[leaveIdx] != theLBbound[leaveIdx])
                       && (enterStat != Desc::P_FREE)
                       && (enterStat != Desc::D_FREE);
      else
         m_numCycle /= 2;

      // setup for updating the copricing vector
      if (coSolveVector2)
         SPxBasis::coSolve(theCoPvec->delta(), *coSolveVector2,
                            unitVecs[leaveIdx], *coSolveVector2rhs);
      else
         SPxBasis::coSolve(theCoPvec->delta(), unitVecs[leaveIdx]);

      (*theCoPrhs)[leaveIdx] = enterRO;
      theCoPvec->value() = (enterRO - enterPric)
                           / theFvec->delta()[leaveIdx];

      if (theCoPvec->value() > epsilon() || theCoPvec->value() < -epsilon())
      {
         if (pricing() == FULL)
         {
            thePvec->value() = theCoPvec->value();
            setupPupdate();
         }
         doPupdate();
      }
      assert(thePvec->isConsistent());
      assert(theCoPvec->isConsistent());

      Real leavebound;             // bound on which leaving variable moves
      getEnterVals2(leaveIdx, enterMax, leavebound);

      //  process entering variable
      theUBbound[leaveIdx] = enterUB;
      theLBbound[leaveIdx] = enterLB;

      //  compute tests:
      updateCoTest();
      if (pricing() == FULL)
         updateTest();


      // update feasibility vectors
      theFvec->value() = leaveVal;
      theFvec->update();
      (*theFvec)[leaveIdx] = enterVal - leaveVal;

      if (leavebound > epsilon() || leavebound < -epsilon())
         theFrhs->multAdd(-leavebound, baseVec(leaveIdx));

      if (enterVal > epsilon() || enterVal < -epsilon())
         theFrhs->multAdd(enterVal, *enterVec);


      //  change basis matrix
      change(leaveIdx, enterId, enterVec, &(theFvec->delta()));
   }


   /*  No leaving vector could be found that would yield a stable pivot step.
    */
   else if (leaveVal != -enterMax)
   {
      rejectEnter(enterId, 0.01*enterTest - 2*delta(), enterStat);
      change(-1, enterId, enterVec);
   }

   /*
       No leaving vector has been selected from the basis. However, if the
       shift amount for |fVec| is bounded, we are in the case, that the
       entering variable is moved from one bound to its other, before any of
       the basis feasibility variables reaches their bound. This may only
       happen in primal/columnwise case with upper and lower bounds on
       variables.
    */
   else if (leaveVal < infinity && leaveVal > -infinity)
   {
      assert(rep() == COLUMN);
      assert(leaveVal == -enterMax);

      change(-1, enterId, enterVec);

      theFvec->value() = leaveVal;
      theFvec->update();

      ungetEnterVal(enterId, enterStat, leaveVal, *enterVec);
   }

   /*
       No variable could be selected to leave the basis and even the entering
       variable is unbounded --- this is a failiour.
    */
   else
   {
      /* The following line originally was in the "lastUpdate() > 1" case;
         we need it in the INFEASIBLE/UNBOUNDED case, too, to have the
         basis descriptor at the correct size. */
      rejectEnter(enterId, enterTest, enterStat);
      if (lastUpdate() > 1)
      {
         std::cout << "AHA!!" << std::endl;
         factorize();
         return enter(enterId);
      }

      SPxId none;
      change(-1, none, 0);
      if (rep() != COLUMN)
         setStatus(SPxBasis::INFEASIBLE);
      else
         setStatus(SPxBasis::UNBOUNDED);
      return 0;
   }

   return 1;
}
} // namespace soplex

//-----------------------------------------------------------------------------
//Emacs Local Variables:
//Emacs mode:c++
//Emacs c-basic-offset:3
//Emacs tab-width:8
//Emacs indent-tabs-mode:nil
//Emacs End:
//-----------------------------------------------------------------------------

invivoenter__ZN6soplex6SoPlex10updateTestEv_204

Real time: 6.912319e+10 -- Predicted time: 6.812312e+10

Callcount: 10000

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5162	### 1	11869.1447732353	1.153320e+05	1.055720e+05	8.4
4325	### 2	6056.62337924111	1.102169e+07	1.129162e+07	2.3

                /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1997-1999 Roland Wunderling                              */
/*                  1997-2002 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.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef SPEC_CPU
#pragma ident "@(#) $Id: ssvector.cpp,v 1.33 2005/01/12 12:00:03 bzfkocht Exp $"
#endif

#include <iostream>
#include <iomanip>
#include <assert.h>

#include "spxdefines.h"
#include "ssvector.h"
#include "svset.h"
#include "message.h"

/**@file ssvector.cpp
 * @todo There is a lot pointer arithmetic done here. It is not clear if
 *       this is an advantage at all. See all the function int() casts.
 * @todo Several operations like maxAbs could setup the vector while
 *       computing there result.
 */
namespace soplex
{
#define MARKER   1e-100

static const Real shortProductFactor = 0.5;

void SSVector::setMax(int newmax)
{
   assert(idx    != 0);
   assert(newmax != 0);
   assert(newmax >= IdxSet::size());

   // len = (newmax < IdxSet::max()) ? IdxSet::max() : newmax;
   len = newmax;

   spx_realloc(idx, len);
}

void SSVector::reDim (int newdim)
{
   for (int i = IdxSet::size() - 1; i >= 0; --i)
      if (index(i) >= newdim)
         remove(i);
   DVector::reDim(newdim);
   setMax(DVector::memSize() + 1);
   assert(isConsistent());
}

void SSVector::reMem(int newsize)
{
   DVector::reSize(newsize);
   assert(isConsistent());
   setMax(DVector::memSize() + 1);
}

void SSVector::clear ()
{
   if (isSetup())
   {
      for(int i = 0; i < num; ++i)
         val[idx[i]] = 0.0;
   }
   else
      Vector::clear();

   IdxSet::clear();
   setupStatus = true;
   assert(isConsistent());
}

void SSVector::setValue(int i, Real x)
{
   assert(i >= 0 && i < DVector::dim());

   if (isSetup())
   {
      int n = number(i);

      if (n < 0)
      {
         if (isNotZero(x, epsilon))
            IdxSet::add(1, &i);
      }
      else if (x == 0)
         clearNum(n);
   }
   val[i] = x;

   assert(isConsistent());
}

#ifdef USE_OLD // old version
void SSVector::setup()
{
   if (!isSetup())
   {
      IdxSet::clear();

      // #define      TWO_LOOPS
#ifdef  TWO_LOOPS
      int i = 0;
      int n = 0;
      int* id = idx;
      Real* v = val;
      const Real* end = val + dim();

      while (v < end)
      {
         id[n] = i++;
         n += (*v++ != 0);
      }

      Real x;
      int* ii = idx;
      int* last = idx + n;
      v = val;

      for (; id < last; ++id)
      {
         x = v[*id];
         if (isNotZero(x, epsilon))
            *ii++ = *id;
         else
            v[*id] = 0;
      }
      num = ii - idx;

#else

      if (dim() <= 1)
      {
         if (dim())
         {
            if (isNotZero(*val, epsilon))
               IdxSet::add(0);
            else
               *val = 0;
         }
      }
      else
      {
         int* ii = idx;
         Real* v = val;
         Real* end = v + dim() - 1;

         /* setze weissen Elefanten */
         Real last = *end;
         *end = MARKER;

         /* erstes element extra */
         if (isNotZero(*v, epsilon))
            *ii++ = 0;
         else
            *v = 0.0;

         for(;;)
         {
            while (*++v == 0.0)
               ;
            if (isNotZero(*v, epsilon))
            {
               *ii++ = int(v - val);
            }
            else
            {
               *v = 0.0;
               if (v == end)
                  break;
            }
         }
         /* fange weissen Elefanten wieder ein */
         if (isNotZero(last, epsilon))
         {
            *v = last;
            *ii++ = dim() - 1;
         }
         else
            *v = 0;

         num = int(ii - idx);
      }

#endif

      setupStatus = true;
      assert(isConsistent());
   }
}
#else // new version, not yet fully testet
void SSVector::setup()
{
   if (!isSetup())
   {
      IdxSet::clear();

      num = 0;

      for(int i = 0; i < dim(); ++i)
      {
         if (val[i] != 0.0)
         {
            if (isZero(val[i], epsilon))
               val[i] = 0.0;
            else
            {
               idx[num] = i;
               num++;
            }
         }
      }
      setupStatus = true;
      assert(isConsistent());
   }
}
#endif

SSVector& SSVector::operator+=(const Vector& vec)
{
   Vector::operator+=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

SSVector& SSVector::operator+=(const SVector& vec)
{
   Vector::operator+=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator+=(const SSVector& vec)
{
   for (int i = vec.size() - 1; i >= 0; --i)
      val[vec.index(i)] += vec.value(i);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#else
SSVector& SSVector::operator+=(const SSVector& vec)
{
   for (int i = 0; i < vec.size(); ++i)
      val[vec.index(i)] += vec.value(i);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#endif

SSVector& SSVector::operator-=(const Vector& vec)
{
   Vector::operator-=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

SSVector& SSVector::operator-=(const SVector& vec)
{
   Vector::operator-=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator-=(const SSVector& vec)
{
   if (vec.isSetup())
   {
      for (int i = vec.size() - 1; i >= 0; --i)
         val[vec.index(i)] -= vec.value(i);
   }
   else
   {
      Vector::operator-=(Vector(vec));
   }

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#else
SSVector& SSVector::operator-=(const SSVector& vec)
{
   if (vec.isSetup())
   {
      for (int i = 0; i < vec.size(); ++i)
         val[vec.index(i)] -= vec.value(i);
   }
   else
   {
      Vector::operator-=(Vector(vec));
   }

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#endif

#ifdef USE_OLD // old version
SSVector& SSVector::operator*=(Real x)
{
   assert(isSetup());

   for (int i = size() - 1; i >= 0; --i)
      val[index(i)] *= x;
   assert(isConsistent());
   return *this;
}
#else
SSVector& SSVector::operator*=(Real x)
{
   assert(isSetup());

   for (int i = 0; i < size(); ++i)
      val[index(i)] *= x;

   assert(isConsistent());

   return *this;
}
#endif

#ifdef USE_OLD // old
Real SSVector::maxAbs() const
{
   if (isSetup())
   {
      int* i = idx;
      int* end = idx + num;
      Real* v = val;
      Real absval = 0.0;

      for (; i < end; ++i)
      {
         Real x = v[*i];
         if (fabs(x) > absval)
            absval = fabs(x);
      }
      return absval;
   }
   else
      return Vector::maxAbs();
}
#else // new, not fully tested
Real SSVector::maxAbs() const
{
   if (isSetup())
   {
      Real maxabs = 0.0;

      for(int i = 0; i < num; ++i)
      {
         Real x = fabs(val[idx[i]]);

         if (x > maxabs)
            maxabs = x;
      }
      return maxabs;
   }
   else
      return Vector::maxAbs();
}
#endif // !0

Real SSVector::length2() const
{
   Real x = 0.0;

   if (isSetup())
   {
      for(int i = 0; i < num; ++i)
         x += val[idx[i]] * val[idx[i]];
   }
   else
      x = Vector::length2();

   return x;
}

Real SSVector::length() const
{
   return sqrt(length2());
}

#if 0 // buggy and not used
/* @todo check if really not used or if the Vector version is used instead.
 */
SSVector& SSVector::multAdd(Real xx, const SSVector& svec)
{
   if (svec.isSetup())
   {
      if (isSetup())
      {
         int i, j;
         Real x;

         for (i = svec.size() - 1; i >= 0; --i)
         {
            j = svec.index(i);
            if (val[j])
            {
               x = val[j] + xx * svec.value(i);
               if (isNotZero(x, epsilon))
                  val[j] = x;
               else
               {
                  val[j] = 0;
                  for (--i; i >= 0; --i)
                     val[svec.index(i)] += xx * svec.value(i);
                  unSetup();
                  break;
               }
            }
            else
            {
               x = xx * svec.value(i);
               if (isNotZero(x, epsilon))
               {
                  val[j] = x;
                  addIdx(j);
               }
            }
         }
      }
      else
         Vector::multAdd(xx, svec);
   }
   else
   {
      /**@todo this code does not work, because in is never something
       *       added to v. Also the idx will not be setup correctly
       *       Fortunately the whole function seems not to be called
       *       at all. 
       */
      assert(false);

      Real y;
      int* ii = idx;
      Real* v = val;
      Real* rv = static_cast<Real*>(svec.val);
      Real* last = rv + svec.dim() - 1;
      Real x = *last;

      *last = MARKER;
      for(;;)
      {
         while (!*rv)
         {
            ++rv;
            ++v;
         }
         y = *rv++ * xx;
         if (isNotZero(y, epsilon))
         {
            *ii++ = int(v - val);
            *v++ = y;
         }
         else if (rv == last)
            break;
         else
            v++;
      }
      *rv = x;

      x *= xx;
      if (isNotZero(x, epsilon))
      {
         *ii++ = int(v - val);
         *v = x;
      }
      num = int(ii - idx);

      setupStatus = true;
   }

   assert(isConsistent());
   return *this;
}
#endif // 0

/* @todo This function does not look good. MARKER is set but never really used.
 */
SSVector& SSVector::multAdd(Real xx, const SVector& svec)
{
   if (isSetup())
   {
      int i, j;
      Real x;
      Real* v = val;
      int adjust = 0;

      for (i = svec.size() - 1; i >= 0; --i)
      {
         j = svec.index(i);
         if (v[j])
         {
            x = v[j] + xx * svec.value(i);
            if (isNotZero(x, epsilon))
               v[j] = x;
            else
            {
               adjust = 1;
               v[j] = MARKER;
            }
         }
         else
         {
            x = xx * svec.value(i);
            if (isNotZero(x, epsilon))
            {
               v[j] = x;
               addIdx(j);
            }
         }
      }

      if (adjust)
      {
         int* iptr = idx;
         int* iiptr = idx;
         int* endptr = idx + num;
         for (; iptr < endptr; ++iptr)
         {
            x = v[*iptr];
            if (isNotZero(x, epsilon))
               *iiptr++ = *iptr;
            else
               v[*iptr] = 0;
         }
         num = int(iiptr - idx);
      }
   }
   else
      Vector::multAdd(xx, svec);

   assert(isConsistent());
   return *this;
}

SSVector& SSVector::multAdd(Real x, const Vector& vec)
{
   Vector::multAdd(x, vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator=(const SSVector& rhs)
{
   assert(rhs.isConsistent());

   if (this != &rhs)
   {
      clear();

      setMax(rhs.max());
      IdxSet::operator=(rhs);
      DVector::reDim(rhs.dim());

      if (rhs.isSetup())
      {
         for (int i = size() - 1; i >= 0; --i)
         {
            int j = index(i);
            val[j] = rhs.val[j];
         }
      }
      else
      {
         int* ii = idx;
         Real* v = val;
         Real* rv = static_cast<Real*>(rhs.val);
         Real* last = rv + rhs.dim() - 1;
         Real x = *last;
         
         *last = MARKER;
         for(;;)
         {
            while (!*rv)
            {
               ++rv;
               ++v;
            }
            if (isNotZero(*rv, epsilon))
            {
               *ii++ = int(v - val);
               *v++ = *rv++;
            }
            else if (rv == last)
               break;
            else
            {
               v++;
               rv++;
            }
         }
         *rv = x;
         
         if (isNotZero(x, epsilon))
         {
            *ii++ = int(v - val);
            *v++ = x;
         }
         num = int(ii - idx);
      }
      setupStatus = true;
   }
   assert(isConsistent());

   return *this;
}
#else // new version
SSVector& SSVector::operator=(const SSVector& rhs)
{
   assert(rhs.isConsistent());

   if (this != &rhs)
   {
      clear();
      epsilon = rhs.epsilon;
      setMax(rhs.max());
      DVector::reDim(rhs.dim());

      if (rhs.isSetup())
      {
         IdxSet::operator=(rhs);

         for(int i = 0; i < size(); ++i)
         {
            int j  = index(i);
            val[j] = rhs.val[j];
         }
      }
      else
      {
         num = 0;

         for(int i = 0; i < rhs.dim(); ++i)
         {
            if (isNotZero(rhs.val[i], epsilon))
            {
               val[i]       = rhs.val[i];
               idx[num]     = i;
               num++;
            }
         }
      }
      setupStatus = true;
   }
   assert(isConsistent());

   return *this;
}
#endif // 0

#ifdef USE_OLD // old version
void SSVector::setup_and_assign(SSVector& rhs)
{
   assert(rhs.isConsistent());

   clear();
   epsilon = rhs.epsilon;
   setMax(rhs.max());
   DVector::reDim(rhs.dim());

   if (rhs.isSetup())
   {
      int i, j;
      IdxSet::operator=(rhs);
      for (i = size() - 1; i >= 0; --i)
      {
         j = index(i);
         val[j] = rhs.val[j];
      }
   }
   else
   {
      int* ri = rhs.idx;
      int* ii = idx;
      Real* rv = rhs.val;
      Real* v = val;
      Real* last = rv + rhs.dim() - 1;
      Real x = *last;

      *last = MARKER;

      for(;;)
      {
         while (!*rv)
         {
            ++rv;
            ++v;
         }
         if (isNotZero(*rv, rhs.epsilon))
         {
            *ri++ = *ii++ = int(v - val);
            *v++ = *rv++;
         }
         else if (rv == last)
            break;
         else
         {
            v++;
            *rv++ = 0;
         }
      }

      if (isNotZero(x, rhs.epsilon))
      {
         *ri++ = *ii++ = int(v - val);
         *v++ = *rv = x;
      }
      else
         *rv = 0;
      num = rhs.num = int(ii - idx);
      rhs.setupStatus = true;
   }
   setupStatus = true;

   assert(isConsistent());
}
#else // new version
/* setup rhs and assign to this
 */
void SSVector::setup_and_assign(SSVector& rhs)
{
   clear();
   epsilon = rhs.epsilon;
   setMax(rhs.max());
   DVector::reDim(rhs.dim());

   if (rhs.isSetup())
   {
      IdxSet::operator=(rhs);
      
      for(int i = 0; i < size(); ++i)
      {
         int j  = index(i);
         val[j] = rhs.val[j];
      }
   }
   else
   {
      num = 0;

      for(int i = 0; i < rhs.dim(); ++i)
      {
         if (rhs.val[i] != 0.0)
         {
            if (isNotZero(rhs.val[i], epsilon))
            {
               rhs.idx[num] = i;
               idx[num]     = i;
               val[i]       = rhs.val[i];
               num++;
            }
            else
            {
               rhs.val[i] = 0.0;
            }
         }
      }
      rhs.num         = num;
      rhs.setupStatus = true;
   }
   setupStatus = true;

   assert(rhs.isConsistent());
   assert(isConsistent());
}
#endif // 0


SSVector& SSVector::operator=(const SVector& rhs)
{
   clear();
   return assign(rhs);
}

#ifdef USE_OLD // old version (buggy or optimization dependend)
SSVector& SSVector::assign(const SVector& rhs)
{
   assert(rhs.dim() <= Vector::dim());

   const SVector::Element* e = rhs.m_elem;
   int* p = idx;
   int i = rhs.size();

   while (i--)
   {
      val[*p = e->idx] = e->val;
      p += ((e++)->val != 0);
   }
   num = int(p - idx);
   setupStatus = true;

   assert(isConsistent());
   return *this;
}

#else // new version not yet fully testet, could be put into operator=()
SSVector& SSVector::assign(const SVector& rhs)
{
   assert(rhs.dim() <= Vector::dim());

   num = 0;

   for(int i = 0; i < rhs.size(); ++i)
   {
      int  k = rhs.index(i);
      Real v = rhs.value(i);

      if (isZero(v, epsilon))
         val[k] = 0.0;
      else
      {
         val[k]     = v;
         idx[num++] = k;
      }
   }
   setupStatus = true;

   assert(isConsistent());

   return *this;
}
#endif

SSVector& SSVector::assign2product1(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   const Real* vl = x.val;
   const int* xi = x.idx;

   int* ii = idx;
   SVector* svec = const_cast<SVector*>( & A[*xi] );
   const SVector::Element* e = &(svec->element(0));
   const SVector::Element* last = e + (num = svec->size());
   Real* v = val;
   Real y = vl[*xi];

   for (; e < last; ++e)
      v[ *ii++ = e->idx ] = y * e->val;

   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productShort(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   int i;
   int j;
   Real y;
   int* ii                      = idx;
   const Real* vl               = x.val;
   const int*  xi               = x.idx;
   SVector* svec                = const_cast<SVector*>( & A[*xi] );
   const SVector::Element* e    = &(svec->element(0));
   const SVector::Element* last = e + (num = svec->size());
   Real* v                      = val;
   Real xx                      = vl[*xi++];

   for (; e < last; ++e)
   {
      y      = xx * e->val;
      *ii    = e->idx;
      v[*ii] = y;
      ii    += (y != 0) ? 1 : 0;
   }

   for (i = x.size(); --i > 0;)
   {
      xx = vl[*xi];
      svec = const_cast<SVector*>( & A[*xi++] );
      e = &(svec->element(0));
      for (int k = svec->size(); --k >= 0;)
      {
         j   = e->idx;
         *ii = j;
         y   = v[j];
         ii  += (y == 0) ? 1 : 0;
         y   += xx * e->val;
         e++;
         v[j] = (y != 0) ? y : MARKER;
      }
   }

   int* is = idx;
   int* it = idx;
   for (; is < ii; ++is)
   {
      if (isNotZero(v[*is], epsilon))
         *it++ = *is;
      else
         v[*is] = 0;
   }
   num = int(it - idx);

   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productFull(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   int i;
   const Real* vl = x.val;
   const int* xi = x.idx;

   SVector* svec;
   const SVector::Element* elem;
   const SVector::Element* last;
   Real y;
   Real* v = val;

   for (i = x.size(); i-- > 0; ++xi)
   {
      svec = const_cast<SVector*>( & A[*xi] );
      elem = &(svec->element(0));
      last = elem + svec->size();
      y = vl[*xi];
      for (; elem < last; ++elem)
         v[elem->idx] += y * elem->val;
   }
   return *this;
}

SSVector& SSVector::assign2product4setup(const SVSet& A, const SSVector& x)
{
   assert(A.num() == x.dim());

   assert(x.isSetup());

   clear();

   if (x.size() == 1)
   {
      assign2product1(A, x);
      setupStatus = true;
   }

   else if (Real(x.size())*A.memSize() <= shortProductFactor*dim()*A.num()
             && isSetup())
   {
      assign2productShort(A, x);
      setupStatus = true;
   }

   else
   {
      assign2productFull(A, x);
      setupStatus = false;
   }
   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2product(const SSVector& x, const SVSet& A)
{
   assert(A.num() == dim());

   Real y;

   clear();

   for (int i = dim(); i-- > 0;)
   {
      y = A[i] * x;

      if (isNotZero(y, epsilon))
      {
         val[i] = y;
         IdxSet::addIdx(i);
      }
   }
   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productAndSetup(const SVSet& A, SSVector& x)
{
   if (x.isSetup())
      return assign2product4setup(A, x);

   SVector* svec;
   const SVector::Element* elem;
   const SVector::Element* last;
   Real y;
   Real* v = val;
   int* xi = x.idx;
   Real* xv = x.val;
   Real* end = xv + x.dim() - 1;

   /* setze weissen Elefanten */
   Real lastval = *end;
   *end = MARKER;

   for(;;)
   {
      while (!*xv)
         ++xv;
      if (isNotZero(*xv, epsilon))
      {
         y = *xv;
         svec = const_cast<SVector*>( & A[ *xi++ = int(xv - x.val) ] );
         elem = &(svec->element(0));
         last = elem + svec->size();
         for (; elem < last; ++elem)
            v[elem->idx] += y * elem->val;
      }
      else
      {
         *xv = 0;
         if (xv == end)
            break;
      }
      xv++;
   }

   /* fange weissen Elefanten wieder ein */
   if (isNotZero(lastval, epsilon))
   {
      y = *xv = lastval;
      svec = const_cast<SVector*>( & A[ *xi++ = int(xv - x.val) ] );
      elem = &(svec->element(0));
      last = elem + svec->size();
      for (; elem < last; ++elem)
         v[elem->idx] += y * elem->val;
   }
   else
      *xv = 0;

   x.num = int(xi - x.idx);
   x.setupStatus = true;
   setupStatus = false;

   assert(isConsistent());

   return *this;
}

bool SSVector::isConsistent() const
{
   if (Vector::dim() > IdxSet::max())
      return MSGinconsistent("SSVector");

   if (Vector::dim() < IdxSet::dim())
      return MSGinconsistent("SSVector");

   if (isSetup())
   {
      for (int i = 0; i < Vector::dim(); ++i)
      {
         int j = number(i);

         if (j < 0 && fabs(val[i]) > 0.0) 
         {
            std::cerr << "i= " << i << " idx= " << j << " val= " << std::setprecision(16) << val[i] << std::endl;
            return MSGinconsistent("SSVector");
         }
      }
   }
   return DVector::isConsistent() && IdxSet::isConsistent();
}
} // namespace soplex

//-----------------------------------------------------------------------------
//Emacs Local Variables:
//Emacs mode:c++
//Emacs c-basic-offset:3
//Emacs tab-width:8
//Emacs indent-tabs-mode:nil
//Emacs End:
//-----------------------------------------------------------------------------

invivossvector__ZN6soplex8SSVector5setupEv_212

Real time: 4.007571e+10 -- Predicted time: 3.799148e+10

Callcount: 6199

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
4201	### 1	3271.80613980497	6.106438e+06	6.260264e+06	2.4
2298	### 2	1426.46504118404	6.128241e+06	6.447644e+06	4.9
279	### 3	383.716682762984	6.413405e+06	7.011956e+06	8.5
5461	### 4	1117.03202241155	6.096285e+06	6.663708e+06	8.5

                /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1997-1999 Roland Wunderling                              */
/*                  1997-2002 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.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef SPEC_CPU
#pragma ident "@(#) $Id: ssvector.cpp,v 1.33 2005/01/12 12:00:03 bzfkocht Exp $"
#endif

#include <iostream>
#include <iomanip>
#include <assert.h>

#include "spxdefines.h"
#include "ssvector.h"
#include "svset.h"
#include "message.h"

/**@file ssvector.cpp
 * @todo There is a lot pointer arithmetic done here. It is not clear if
 *       this is an advantage at all. See all the function int() casts.
 * @todo Several operations like maxAbs could setup the vector while
 *       computing there result.
 */
namespace soplex
{
#define MARKER   1e-100

static const Real shortProductFactor = 0.5;

void SSVector::setMax(int newmax)
{
   assert(idx    != 0);
   assert(newmax != 0);
   assert(newmax >= IdxSet::size());

   // len = (newmax < IdxSet::max()) ? IdxSet::max() : newmax;
   len = newmax;

   spx_realloc(idx, len);
}

void SSVector::reDim (int newdim)
{
   for (int i = IdxSet::size() - 1; i >= 0; --i)
      if (index(i) >= newdim)
         remove(i);
   DVector::reDim(newdim);
   setMax(DVector::memSize() + 1);
   assert(isConsistent());
}

void SSVector::reMem(int newsize)
{
   DVector::reSize(newsize);
   assert(isConsistent());
   setMax(DVector::memSize() + 1);
}

void SSVector::clear ()
{
   if (isSetup())
   {
      for(int i = 0; i < num; ++i)
         val[idx[i]] = 0.0;
   }
   else
      Vector::clear();

   IdxSet::clear();
   setupStatus = true;
   assert(isConsistent());
}

void SSVector::setValue(int i, Real x)
{
   assert(i >= 0 && i < DVector::dim());

   if (isSetup())
   {
      int n = number(i);

      if (n < 0)
      {
         if (isNotZero(x, epsilon))
            IdxSet::add(1, &i);
      }
      else if (x == 0)
         clearNum(n);
   }
   val[i] = x;

   assert(isConsistent());
}

#ifdef USE_OLD // old version
void SSVector::setup()
{
   if (!isSetup())
   {
      IdxSet::clear();

      // #define      TWO_LOOPS
#ifdef  TWO_LOOPS
      int i = 0;
      int n = 0;
      int* id = idx;
      Real* v = val;
      const Real* end = val + dim();

      while (v < end)
      {
         id[n] = i++;
         n += (*v++ != 0);
      }

      Real x;
      int* ii = idx;
      int* last = idx + n;
      v = val;

      for (; id < last; ++id)
      {
         x = v[*id];
         if (isNotZero(x, epsilon))
            *ii++ = *id;
         else
            v[*id] = 0;
      }
      num = ii - idx;

#else

      if (dim() <= 1)
      {
         if (dim())
         {
            if (isNotZero(*val, epsilon))
               IdxSet::add(0);
            else
               *val = 0;
         }
      }
      else
      {
         int* ii = idx;
         Real* v = val;
         Real* end = v + dim() - 1;

         /* setze weissen Elefanten */
         Real last = *end;
         *end = MARKER;

         /* erstes element extra */
         if (isNotZero(*v, epsilon))
            *ii++ = 0;
         else
            *v = 0.0;

         for(;;)
         {
            while (*++v == 0.0)
               ;
            if (isNotZero(*v, epsilon))
            {
               *ii++ = int(v - val);
            }
            else
            {
               *v = 0.0;
               if (v == end)
                  break;
            }
         }
         /* fange weissen Elefanten wieder ein */
         if (isNotZero(last, epsilon))
         {
            *v = last;
            *ii++ = dim() - 1;
         }
         else
            *v = 0;

         num = int(ii - idx);
      }

#endif

      setupStatus = true;
      assert(isConsistent());
   }
}
#else // new version, not yet fully testet
void SSVector::setup()
{
   if (!isSetup())
   {
      IdxSet::clear();

      num = 0;

      for(int i = 0; i < dim(); ++i)
      {
         if (val[i] != 0.0)
         {
            if (isZero(val[i], epsilon))
               val[i] = 0.0;
            else
            {
               idx[num] = i;
               num++;
            }
         }
      }
      setupStatus = true;
      assert(isConsistent());
   }
}
#endif

SSVector& SSVector::operator+=(const Vector& vec)
{
   Vector::operator+=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

SSVector& SSVector::operator+=(const SVector& vec)
{
   Vector::operator+=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator+=(const SSVector& vec)
{
   for (int i = vec.size() - 1; i >= 0; --i)
      val[vec.index(i)] += vec.value(i);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#else
SSVector& SSVector::operator+=(const SSVector& vec)
{
   for (int i = 0; i < vec.size(); ++i)
      val[vec.index(i)] += vec.value(i);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#endif

SSVector& SSVector::operator-=(const Vector& vec)
{
   Vector::operator-=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

SSVector& SSVector::operator-=(const SVector& vec)
{
   Vector::operator-=(vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator-=(const SSVector& vec)
{
   if (vec.isSetup())
   {
      for (int i = vec.size() - 1; i >= 0; --i)
         val[vec.index(i)] -= vec.value(i);
   }
   else
   {
      Vector::operator-=(Vector(vec));
   }

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#else
SSVector& SSVector::operator-=(const SSVector& vec)
{
   if (vec.isSetup())
   {
      for (int i = 0; i < vec.size(); ++i)
         val[vec.index(i)] -= vec.value(i);
   }
   else
   {
      Vector::operator-=(Vector(vec));
   }

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}
#endif

#ifdef USE_OLD // old version
SSVector& SSVector::operator*=(Real x)
{
   assert(isSetup());

   for (int i = size() - 1; i >= 0; --i)
      val[index(i)] *= x;
   assert(isConsistent());
   return *this;
}
#else
SSVector& SSVector::operator*=(Real x)
{
   assert(isSetup());

   for (int i = 0; i < size(); ++i)
      val[index(i)] *= x;

   assert(isConsistent());

   return *this;
}
#endif

#ifdef USE_OLD // old
Real SSVector::maxAbs() const
{
   if (isSetup())
   {
      int* i = idx;
      int* end = idx + num;
      Real* v = val;
      Real absval = 0.0;

      for (; i < end; ++i)
      {
         Real x = v[*i];
         if (fabs(x) > absval)
            absval = fabs(x);
      }
      return absval;
   }
   else
      return Vector::maxAbs();
}
#else // new, not fully tested
Real SSVector::maxAbs() const
{
   if (isSetup())
   {
      Real maxabs = 0.0;

      for(int i = 0; i < num; ++i)
      {
         Real x = fabs(val[idx[i]]);

         if (x > maxabs)
            maxabs = x;
      }
      return maxabs;
   }
   else
      return Vector::maxAbs();
}
#endif // !0

Real SSVector::length2() const
{
   Real x = 0.0;

   if (isSetup())
   {
      for(int i = 0; i < num; ++i)
         x += val[idx[i]] * val[idx[i]];
   }
   else
      x = Vector::length2();

   return x;
}

Real SSVector::length() const
{
   return sqrt(length2());
}

#if 0 // buggy and not used
/* @todo check if really not used or if the Vector version is used instead.
 */
SSVector& SSVector::multAdd(Real xx, const SSVector& svec)
{
   if (svec.isSetup())
   {
      if (isSetup())
      {
         int i, j;
         Real x;

         for (i = svec.size() - 1; i >= 0; --i)
         {
            j = svec.index(i);
            if (val[j])
            {
               x = val[j] + xx * svec.value(i);
               if (isNotZero(x, epsilon))
                  val[j] = x;
               else
               {
                  val[j] = 0;
                  for (--i; i >= 0; --i)
                     val[svec.index(i)] += xx * svec.value(i);
                  unSetup();
                  break;
               }
            }
            else
            {
               x = xx * svec.value(i);
               if (isNotZero(x, epsilon))
               {
                  val[j] = x;
                  addIdx(j);
               }
            }
         }
      }
      else
         Vector::multAdd(xx, svec);
   }
   else
   {
      /**@todo this code does not work, because in is never something
       *       added to v. Also the idx will not be setup correctly
       *       Fortunately the whole function seems not to be called
       *       at all. 
       */
      assert(false);

      Real y;
      int* ii = idx;
      Real* v = val;
      Real* rv = static_cast<Real*>(svec.val);
      Real* last = rv + svec.dim() - 1;
      Real x = *last;

      *last = MARKER;
      for(;;)
      {
         while (!*rv)
         {
            ++rv;
            ++v;
         }
         y = *rv++ * xx;
         if (isNotZero(y, epsilon))
         {
            *ii++ = int(v - val);
            *v++ = y;
         }
         else if (rv == last)
            break;
         else
            v++;
      }
      *rv = x;

      x *= xx;
      if (isNotZero(x, epsilon))
      {
         *ii++ = int(v - val);
         *v = x;
      }
      num = int(ii - idx);

      setupStatus = true;
   }

   assert(isConsistent());
   return *this;
}
#endif // 0

/* @todo This function does not look good. MARKER is set but never really used.
 */
SSVector& SSVector::multAdd(Real xx, const SVector& svec)
{
   if (isSetup())
   {
      int i, j;
      Real x;
      Real* v = val;
      int adjust = 0;

      for (i = svec.size() - 1; i >= 0; --i)
      {
         j = svec.index(i);
         if (v[j])
         {
            x = v[j] + xx * svec.value(i);
            if (isNotZero(x, epsilon))
               v[j] = x;
            else
            {
               adjust = 1;
               v[j] = MARKER;
            }
         }
         else
         {
            x = xx * svec.value(i);
            if (isNotZero(x, epsilon))
            {
               v[j] = x;
               addIdx(j);
            }
         }
      }

      if (adjust)
      {
         int* iptr = idx;
         int* iiptr = idx;
         int* endptr = idx + num;
         for (; iptr < endptr; ++iptr)
         {
            x = v[*iptr];
            if (isNotZero(x, epsilon))
               *iiptr++ = *iptr;
            else
               v[*iptr] = 0;
         }
         num = int(iiptr - idx);
      }
   }
   else
      Vector::multAdd(xx, svec);

   assert(isConsistent());
   return *this;
}

SSVector& SSVector::multAdd(Real x, const Vector& vec)
{
   Vector::multAdd(x, vec);

   if (isSetup())
   {
      setupStatus = false;
      setup();
   }
   return *this;
}

#ifdef USE_OLD // old version
SSVector& SSVector::operator=(const SSVector& rhs)
{
   assert(rhs.isConsistent());

   if (this != &rhs)
   {
      clear();

      setMax(rhs.max());
      IdxSet::operator=(rhs);
      DVector::reDim(rhs.dim());

      if (rhs.isSetup())
      {
         for (int i = size() - 1; i >= 0; --i)
         {
            int j = index(i);
            val[j] = rhs.val[j];
         }
      }
      else
      {
         int* ii = idx;
         Real* v = val;
         Real* rv = static_cast<Real*>(rhs.val);
         Real* last = rv + rhs.dim() - 1;
         Real x = *last;
         
         *last = MARKER;
         for(;;)
         {
            while (!*rv)
            {
               ++rv;
               ++v;
            }
            if (isNotZero(*rv, epsilon))
            {
               *ii++ = int(v - val);
               *v++ = *rv++;
            }
            else if (rv == last)
               break;
            else
            {
               v++;
               rv++;
            }
         }
         *rv = x;
         
         if (isNotZero(x, epsilon))
         {
            *ii++ = int(v - val);
            *v++ = x;
         }
         num = int(ii - idx);
      }
      setupStatus = true;
   }
   assert(isConsistent());

   return *this;
}
#else // new version
SSVector& SSVector::operator=(const SSVector& rhs)
{
   assert(rhs.isConsistent());

   if (this != &rhs)
   {
      clear();
      epsilon = rhs.epsilon;
      setMax(rhs.max());
      DVector::reDim(rhs.dim());

      if (rhs.isSetup())
      {
         IdxSet::operator=(rhs);

         for(int i = 0; i < size(); ++i)
         {
            int j  = index(i);
            val[j] = rhs.val[j];
         }
      }
      else
      {
         num = 0;

         for(int i = 0; i < rhs.dim(); ++i)
         {
            if (isNotZero(rhs.val[i], epsilon))
            {
               val[i]       = rhs.val[i];
               idx[num]     = i;
               num++;
            }
         }
      }
      setupStatus = true;
   }
   assert(isConsistent());

   return *this;
}
#endif // 0

#ifdef USE_OLD // old version
void SSVector::setup_and_assign(SSVector& rhs)
{
   assert(rhs.isConsistent());

   clear();
   epsilon = rhs.epsilon;
   setMax(rhs.max());
   DVector::reDim(rhs.dim());

   if (rhs.isSetup())
   {
      int i, j;
      IdxSet::operator=(rhs);
      for (i = size() - 1; i >= 0; --i)
      {
         j = index(i);
         val[j] = rhs.val[j];
      }
   }
   else
   {
      int* ri = rhs.idx;
      int* ii = idx;
      Real* rv = rhs.val;
      Real* v = val;
      Real* last = rv + rhs.dim() - 1;
      Real x = *last;

      *last = MARKER;

      for(;;)
      {
         while (!*rv)
         {
            ++rv;
            ++v;
         }
         if (isNotZero(*rv, rhs.epsilon))
         {
            *ri++ = *ii++ = int(v - val);
            *v++ = *rv++;
         }
         else if (rv == last)
            break;
         else
         {
            v++;
            *rv++ = 0;
         }
      }

      if (isNotZero(x, rhs.epsilon))
      {
         *ri++ = *ii++ = int(v - val);
         *v++ = *rv = x;
      }
      else
         *rv = 0;
      num = rhs.num = int(ii - idx);
      rhs.setupStatus = true;
   }
   setupStatus = true;

   assert(isConsistent());
}
#else // new version
/* setup rhs and assign to this
 */
void SSVector::setup_and_assign(SSVector& rhs)
{
   clear();
   epsilon = rhs.epsilon;
   setMax(rhs.max());
   DVector::reDim(rhs.dim());

   if (rhs.isSetup())
   {
      IdxSet::operator=(rhs);
      
      for(int i = 0; i < size(); ++i)
      {
         int j  = index(i);
         val[j] = rhs.val[j];
      }
   }
   else
   {
      num = 0;

      for(int i = 0; i < rhs.dim(); ++i)
      {
         if (rhs.val[i] != 0.0)
         {
            if (isNotZero(rhs.val[i], epsilon))
            {
               rhs.idx[num] = i;
               idx[num]     = i;
               val[i]       = rhs.val[i];
               num++;
            }
            else
            {
               rhs.val[i] = 0.0;
            }
         }
      }
      rhs.num         = num;
      rhs.setupStatus = true;
   }
   setupStatus = true;

   assert(rhs.isConsistent());
   assert(isConsistent());
}
#endif // 0


SSVector& SSVector::operator=(const SVector& rhs)
{
   clear();
   return assign(rhs);
}

#ifdef USE_OLD // old version (buggy or optimization dependend)
SSVector& SSVector::assign(const SVector& rhs)
{
   assert(rhs.dim() <= Vector::dim());

   const SVector::Element* e = rhs.m_elem;
   int* p = idx;
   int i = rhs.size();

   while (i--)
   {
      val[*p = e->idx] = e->val;
      p += ((e++)->val != 0);
   }
   num = int(p - idx);
   setupStatus = true;

   assert(isConsistent());
   return *this;
}

#else // new version not yet fully testet, could be put into operator=()
SSVector& SSVector::assign(const SVector& rhs)
{
   assert(rhs.dim() <= Vector::dim());

   num = 0;

   for(int i = 0; i < rhs.size(); ++i)
   {
      int  k = rhs.index(i);
      Real v = rhs.value(i);

      if (isZero(v, epsilon))
         val[k] = 0.0;
      else
      {
         val[k]     = v;
         idx[num++] = k;
      }
   }
   setupStatus = true;

   assert(isConsistent());

   return *this;
}
#endif

SSVector& SSVector::assign2product1(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   const Real* vl = x.val;
   const int* xi = x.idx;

   int* ii = idx;
   SVector* svec = const_cast<SVector*>( & A[*xi] );
   const SVector::Element* e = &(svec->element(0));
   const SVector::Element* last = e + (num = svec->size());
   Real* v = val;
   Real y = vl[*xi];

   for (; e < last; ++e)
      v[ *ii++ = e->idx ] = y * e->val;

   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productShort(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   int i;
   int j;
   Real y;
   int* ii                      = idx;
   const Real* vl               = x.val;
   const int*  xi               = x.idx;
   SVector* svec                = const_cast<SVector*>( & A[*xi] );
   const SVector::Element* e    = &(svec->element(0));
   const SVector::Element* last = e + (num = svec->size());
   Real* v                      = val;
   Real xx                      = vl[*xi++];

   for (; e < last; ++e)
   {
      y      = xx * e->val;
      *ii    = e->idx;
      v[*ii] = y;
      ii    += (y != 0) ? 1 : 0;
   }

   for (i = x.size(); --i > 0;)
   {
      xx = vl[*xi];
      svec = const_cast<SVector*>( & A[*xi++] );
      e = &(svec->element(0));
      for (int k = svec->size(); --k >= 0;)
      {
         j   = e->idx;
         *ii = j;
         y   = v[j];
         ii  += (y == 0) ? 1 : 0;
         y   += xx * e->val;
         e++;
         v[j] = (y != 0) ? y : MARKER;
      }
   }

   int* is = idx;
   int* it = idx;
   for (; is < ii; ++is)
   {
      if (isNotZero(v[*is], epsilon))
         *it++ = *is;
      else
         v[*is] = 0;
   }
   num = int(it - idx);

   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productFull(const SVSet& A, const SSVector& x)
{
   assert(x.isSetup());

   int i;
   const Real* vl = x.val;
   const int* xi = x.idx;

   SVector* svec;
   const SVector::Element* elem;
   const SVector::Element* last;
   Real y;
   Real* v = val;

   for (i = x.size(); i-- > 0; ++xi)
   {
      svec = const_cast<SVector*>( & A[*xi] );
      elem = &(svec->element(0));
      last = elem + svec->size();
      y = vl[*xi];
      for (; elem < last; ++elem)
         v[elem->idx] += y * elem->val;
   }
   return *this;
}

SSVector& SSVector::assign2product4setup(const SVSet& A, const SSVector& x)
{
   assert(A.num() == x.dim());

   assert(x.isSetup());

   clear();

   if (x.size() == 1)
   {
      assign2product1(A, x);
      setupStatus = true;
   }

   else if (Real(x.size())*A.memSize() <= shortProductFactor*dim()*A.num()
             && isSetup())
   {
      assign2productShort(A, x);
      setupStatus = true;
   }

   else
   {
      assign2productFull(A, x);
      setupStatus = false;
   }
   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2product(const SSVector& x, const SVSet& A)
{
   assert(A.num() == dim());

   Real y;

   clear();

   for (int i = dim(); i-- > 0;)
   {
      y = A[i] * x;

      if (isNotZero(y, epsilon))
      {
         val[i] = y;
         IdxSet::addIdx(i);
      }
   }
   assert(isConsistent());

   return *this;
}

SSVector& SSVector::assign2productAndSetup(const SVSet& A, SSVector& x)
{
   if (x.isSetup())
      return assign2product4setup(A, x);

   SVector* svec;
   const SVector::Element* elem;
   const SVector::Element* last;
   Real y;
   Real* v = val;
   int* xi = x.idx;
   Real* xv = x.val;
   Real* end = xv + x.dim() - 1;

   /* setze weissen Elefanten */
   Real lastval = *end;
   *end = MARKER;

   for(;;)
   {
      while (!*xv)
         ++xv;
      if (isNotZero(*xv, epsilon))
      {
         y = *xv;
         svec = const_cast<SVector*>( & A[ *xi++ = int(xv - x.val) ] );
         elem = &(svec->element(0));
         last = elem + svec->size();
         for (; elem < last; ++elem)
            v[elem->idx] += y * elem->val;
      }
      else
      {
         *xv = 0;
         if (xv == end)
            break;
      }
      xv++;
   }

   /* fange weissen Elefanten wieder ein */
   if (isNotZero(lastval, epsilon))
   {
      y = *xv = lastval;
      svec = const_cast<SVector*>( & A[ *xi++ = int(xv - x.val) ] );
      elem = &(svec->element(0));
      last = elem + svec->size();
      for (; elem < last; ++elem)
         v[elem->idx] += y * elem->val;
   }
   else
      *xv = 0;

   x.num = int(xi - x.idx);
   x.setupStatus = true;
   setupStatus = false;

   assert(isConsistent());

   return *this;
}

bool SSVector::isConsistent() const
{
   if (Vector::dim() > IdxSet::max())
      return MSGinconsistent("SSVector");

   if (Vector::dim() < IdxSet::dim())
      return MSGinconsistent("SSVector");

   if (isSetup())
   {
      for (int i = 0; i < Vector::dim(); ++i)
      {
         int j = number(i);

         if (j < 0 && fabs(val[i]) > 0.0) 
         {
            std::cerr << "i= " << i << " idx= " << j << " val= " << std::setprecision(16) << val[i] << std::endl;
            return MSGinconsistent("SSVector");
         }
      }
   }
   return DVector::isConsistent() && IdxSet::isConsistent();
}
} // namespace soplex

//-----------------------------------------------------------------------------
//Emacs Local Variables:
//Emacs mode:c++
//Emacs c-basic-offset:3
//Emacs tab-width:8
//Emacs indent-tabs-mode:nil
//Emacs End:
//-----------------------------------------------------------------------------

invivossvector__ZN6soplex8SSVector18assign2productFullERKNS_5SVSetERKS0__977

Real time: 1.273837e+11 -- Predicted time: 1.346080e+11

Callcount: 6199

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
540	### 1	721.731340654192	1.737297e+07	1.687015e+07	2.8
5859	### 2	5411.56186298866	2.255714e+07	2.183124e+07	3.2

450-soplex

Exec Time : 8.433491e+11 cycles

Matching : 39.9%

__invivo__spxsteeppr__ZN6soplex10SPxSteepPR9entered4XENS_5SPxIdEiiiii_399

Real time: 3.997608e+11 -- Predicted time: 7.590169e+09

Callcount: 10000

__invivo__spxsteeppr__ZN6soplex10SPxSteepPR12selectEnterXERdiiii_449

Real time: 1.015031e+11 -- Predicted time: 1.070856e+11

Callcount: 10000

__invivo__spxsteeppr___svector__ZNK6soplex7SVectormlERKNS_6VectorE_293

Real time: 3.189254e+11 -- Predicted time: 1.094412e+11

Callcount: 696945308

__invivo__ssvector__ZN6soplex8SSVector5clearEv_74

Real time: 1.639992e+10 -- Predicted time: 1.690773e+10

Callcount: 90023

__invivo__enter__ZN6soplex6SoPlex10updateTestEv_204

Real time: 6.912319e+10 -- Predicted time: 6.812312e+10

Callcount: 10000

__invivo__ssvector__ZN6soplex8SSVector5setupEv_212

Real time: 4.007571e+10 -- Predicted time: 3.799148e+10

Callcount: 6199

__invivo__ssvector__ZN6soplex8SSVector18assign2productFullERKNS_5SVSetERKS0__977

Real time: 1.273837e+11 -- Predicted time: 1.346080e+11

Callcount: 6199

invivospxsteeppr__ZN6soplex10SPxSteepPR9entered4XENS_5SPxIdEiiiii_399

invivospxsteeppr__ZN6soplex10SPxSteepPR12selectEnterXERdiiii_449

invivospxsteeppr_svectorZNK6soplex7SVectormlERKNS_6VectorE_293

invivossvector__ZN6soplex8SSVector5clearEv_74

invivoenter__ZN6soplex6SoPlex10updateTestEv_204

invivossvector__ZN6soplex8SSVector5setupEv_212

invivossvector__ZN6soplex8SSVector18assign2productFullERKNS_5SVSetERKS0__977