Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpCholeskyBase.cpp

    r2271 r2385  
    2121/*----------------------------------------------------------------------------*/
    2222
    23 
    2423#include "CoinPragma.hpp"
    2524
     
    4241// Default Constructor
    4342//-------------------------------------------------------------------
    44 ClpCholeskyBase::ClpCholeskyBase (int denseThreshold) :
    45      type_(0),
    46      doKKT_(false),
    47      goDense_(0.7),
    48      choleskyCondition_(0.0),
    49      model_(NULL),
    50      numberTrials_(),
    51      numberRows_(0),
    52      status_(0),
    53      rowsDropped_(NULL),
    54      permuteInverse_(NULL),
    55      permute_(NULL),
    56      numberRowsDropped_(0),
    57      sparseFactor_(NULL),
    58      choleskyStart_(NULL),
    59      choleskyRow_(NULL),
    60      indexStart_(NULL),
    61      diagonal_(NULL),
    62      workDouble_(NULL),
    63      link_(NULL),
    64      workInteger_(NULL),
    65      clique_(NULL),
    66      sizeFactor_(0),
    67      sizeIndex_(0),
    68      firstDense_(0),
    69      rowCopy_(NULL),
    70      whichDense_(NULL),
    71      denseColumn_(NULL),
    72      dense_(NULL),
    73     denseThreshold_(denseThreshold)
     43ClpCholeskyBase::ClpCholeskyBase(int denseThreshold)
     44  : type_(0)
     45  , doKKT_(false)
     46  , goDense_(0.7)
     47  , choleskyCondition_(0.0)
     48  , model_(NULL)
     49  , numberTrials_()
     50  , numberRows_(0)
     51  , status_(0)
     52  , rowsDropped_(NULL)
     53  , permuteInverse_(NULL)
     54  , permute_(NULL)
     55  , numberRowsDropped_(0)
     56  , sparseFactor_(NULL)
     57  , choleskyStart_(NULL)
     58  , choleskyRow_(NULL)
     59  , indexStart_(NULL)
     60  , diagonal_(NULL)
     61  , workDouble_(NULL)
     62  , link_(NULL)
     63  , workInteger_(NULL)
     64  , clique_(NULL)
     65  , sizeFactor_(0)
     66  , sizeIndex_(0)
     67  , firstDense_(0)
     68  , rowCopy_(NULL)
     69  , whichDense_(NULL)
     70  , denseColumn_(NULL)
     71  , dense_(NULL)
     72  , denseThreshold_(denseThreshold)
    7473{
    75      memset(integerParameters_, 0, 64 * sizeof(int));
    76      memset(doubleParameters_, 0, 64 * sizeof(double));
     74  memset(integerParameters_, 0, 64 * sizeof(int));
     75  memset(doubleParameters_, 0, 64 * sizeof(double));
    7776}
    7877
     
    8079// Copy constructor
    8180//-------------------------------------------------------------------
    82 ClpCholeskyBase::ClpCholeskyBase (const ClpCholeskyBase & rhs) :
    83      type_(rhs.type_),
    84      doKKT_(rhs.doKKT_),
    85      goDense_(rhs.goDense_),
    86      choleskyCondition_(rhs.choleskyCondition_),
    87      model_(rhs.model_),
    88      numberTrials_(rhs.numberTrials_),
    89      numberRows_(rhs.numberRows_),
    90      status_(rhs.status_),
    91     numberRowsDropped_(rhs.numberRowsDropped_)
     81ClpCholeskyBase::ClpCholeskyBase(const ClpCholeskyBase &rhs)
     82  : type_(rhs.type_)
     83  , doKKT_(rhs.doKKT_)
     84  , goDense_(rhs.goDense_)
     85  , choleskyCondition_(rhs.choleskyCondition_)
     86  , model_(rhs.model_)
     87  , numberTrials_(rhs.numberTrials_)
     88  , numberRows_(rhs.numberRows_)
     89  , status_(rhs.status_)
     90  , numberRowsDropped_(rhs.numberRowsDropped_)
    9291{
    93      rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
    94      permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
    95      permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
    96      sizeFactor_ = rhs.sizeFactor_;
    97      sizeIndex_ = rhs.sizeIndex_;
    98      firstDense_ = rhs.firstDense_;
    99      sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
    100      choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
    101      indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
    102      choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
    103      diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
    104 #if CLP_LONG_CHOLESKY!=1
    105      workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
     92  rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
     93  permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
     94  permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
     95  sizeFactor_ = rhs.sizeFactor_;
     96  sizeIndex_ = rhs.sizeIndex_;
     97  firstDense_ = rhs.firstDense_;
     98  sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
     99  choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
     100  indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
     101  choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
     102  diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
     103#if CLP_LONG_CHOLESKY != 1
     104  workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
    106105#else
    107      // actually long double
    108      workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_), numberRows_));
    109 #endif
    110      link_ = ClpCopyOfArray(rhs.link_, numberRows_);
    111      workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
    112      clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
    113      CoinMemcpyN(rhs.integerParameters_, 64, integerParameters_);
    114      CoinMemcpyN(rhs.doubleParameters_, 64, doubleParameters_);
    115      rowCopy_ = rhs.rowCopy_->clone();
    116      whichDense_ = NULL;
    117      denseColumn_ = NULL;
    118      dense_ = NULL;
    119      denseThreshold_ = rhs.denseThreshold_;
     106  // actually long double
     107  workDouble_ = reinterpret_cast< double * >(ClpCopyOfArray(reinterpret_cast< CoinWorkDouble * >(rhs.workDouble_), numberRows_));
     108#endif
     109  link_ = ClpCopyOfArray(rhs.link_, numberRows_);
     110  workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
     111  clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
     112  CoinMemcpyN(rhs.integerParameters_, 64, integerParameters_);
     113  CoinMemcpyN(rhs.doubleParameters_, 64, doubleParameters_);
     114  rowCopy_ = rhs.rowCopy_->clone();
     115  whichDense_ = NULL;
     116  denseColumn_ = NULL;
     117  dense_ = NULL;
     118  denseThreshold_ = rhs.denseThreshold_;
    120119}
    121120
     
    123122// Destructor
    124123//-------------------------------------------------------------------
    125 ClpCholeskyBase::~ClpCholeskyBase ()
     124ClpCholeskyBase::~ClpCholeskyBase()
    126125{
    127      delete [] rowsDropped_;
    128      delete [] permuteInverse_;
    129      delete [] permute_;
    130      delete [] sparseFactor_;
    131      delete [] choleskyStart_;
    132      delete [] choleskyRow_;
    133      delete [] indexStart_;
    134      delete [] diagonal_;
    135      delete [] workDouble_;
    136      delete [] link_;
    137      delete [] workInteger_;
    138      delete [] clique_;
    139      delete rowCopy_;
    140      delete [] whichDense_;
    141      delete [] denseColumn_;
    142      delete dense_;
     126  delete[] rowsDropped_;
     127  delete[] permuteInverse_;
     128  delete[] permute_;
     129  delete[] sparseFactor_;
     130  delete[] choleskyStart_;
     131  delete[] choleskyRow_;
     132  delete[] indexStart_;
     133  delete[] diagonal_;
     134  delete[] workDouble_;
     135  delete[] link_;
     136  delete[] workInteger_;
     137  delete[] clique_;
     138  delete rowCopy_;
     139  delete[] whichDense_;
     140  delete[] denseColumn_;
     141  delete dense_;
    143142}
    144143
     
    147146//-------------------------------------------------------------------
    148147ClpCholeskyBase &
    149 ClpCholeskyBase::operator=(const ClpCholeskyBase& rhs)
     148ClpCholeskyBase::operator=(const ClpCholeskyBase &rhs)
    150149{
    151      if (this != &rhs) {
    152           type_ = rhs.type_;
    153           doKKT_ = rhs.doKKT_;
    154           goDense_ = rhs.goDense_;
    155           choleskyCondition_ = rhs.choleskyCondition_;
    156           model_ = rhs.model_;
    157           numberTrials_ = rhs.numberTrials_;
    158           numberRows_ = rhs.numberRows_;
    159           status_ = rhs.status_;
    160           numberRowsDropped_ = rhs.numberRowsDropped_;
    161           delete [] rowsDropped_;
    162           delete [] permuteInverse_;
    163           delete [] permute_;
    164           delete [] sparseFactor_;
    165           delete [] choleskyStart_;
    166           delete [] choleskyRow_;
    167           delete [] indexStart_;
    168           delete [] diagonal_;
    169           delete [] workDouble_;
    170           delete [] link_;
    171           delete [] workInteger_;
    172           delete [] clique_;
    173           delete rowCopy_;
    174           delete [] whichDense_;
    175           delete [] denseColumn_;
    176           delete dense_;
    177           rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
    178           permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
    179           permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
    180           sizeFactor_ = rhs.sizeFactor_;
    181           sizeIndex_ = rhs.sizeIndex_;
    182           firstDense_ = rhs.firstDense_;
    183           sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
    184           choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
    185           choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, rhs.sizeFactor_);
    186           indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
    187           choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
    188           diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
    189 #if CLP_LONG_CHOLESKY!=1
    190           workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
     150  if (this != &rhs) {
     151    type_ = rhs.type_;
     152    doKKT_ = rhs.doKKT_;
     153    goDense_ = rhs.goDense_;
     154    choleskyCondition_ = rhs.choleskyCondition_;
     155    model_ = rhs.model_;
     156    numberTrials_ = rhs.numberTrials_;
     157    numberRows_ = rhs.numberRows_;
     158    status_ = rhs.status_;
     159    numberRowsDropped_ = rhs.numberRowsDropped_;
     160    delete[] rowsDropped_;
     161    delete[] permuteInverse_;
     162    delete[] permute_;
     163    delete[] sparseFactor_;
     164    delete[] choleskyStart_;
     165    delete[] choleskyRow_;
     166    delete[] indexStart_;
     167    delete[] diagonal_;
     168    delete[] workDouble_;
     169    delete[] link_;
     170    delete[] workInteger_;
     171    delete[] clique_;
     172    delete rowCopy_;
     173    delete[] whichDense_;
     174    delete[] denseColumn_;
     175    delete dense_;
     176    rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
     177    permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
     178    permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
     179    sizeFactor_ = rhs.sizeFactor_;
     180    sizeIndex_ = rhs.sizeIndex_;
     181    firstDense_ = rhs.firstDense_;
     182    sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
     183    choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
     184    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, rhs.sizeFactor_);
     185    indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
     186    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
     187    diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
     188#if CLP_LONG_CHOLESKY != 1
     189    workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
    191190#else
    192           // actually long double
    193           workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_), numberRows_));
    194 #endif
    195           link_ = ClpCopyOfArray(rhs.link_, numberRows_);
    196           workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
    197           clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
    198           rowCopy_ = rhs.rowCopy_->clone();
    199           whichDense_ = NULL;
    200           denseColumn_ = NULL;
    201           dense_ = NULL;
    202           denseThreshold_ = rhs.denseThreshold_;
    203      }
    204      return *this;
     191    // actually long double
     192    workDouble_ = reinterpret_cast< double * >(ClpCopyOfArray(reinterpret_cast< CoinWorkDouble * >(rhs.workDouble_), numberRows_));
     193#endif
     194    link_ = ClpCopyOfArray(rhs.link_, numberRows_);
     195    workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
     196    clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
     197    rowCopy_ = rhs.rowCopy_->clone();
     198    whichDense_ = NULL;
     199    denseColumn_ = NULL;
     200    dense_ = NULL;
     201    denseThreshold_ = rhs.denseThreshold_;
     202  }
     203  return *this;
    205204}
    206205// reset numberRowsDropped and rowsDropped.
    207 void
    208 ClpCholeskyBase::resetRowsDropped()
     206void ClpCholeskyBase::resetRowsDropped()
    209207{
    210      numberRowsDropped_ = 0;
    211      memset(rowsDropped_, 0, numberRows_);
     208  numberRowsDropped_ = 0;
     209  memset(rowsDropped_, 0, numberRows_);
    212210}
    213211/* Uses factorization to solve. - given as if KKT.
    214212   region1 is rows+columns, region2 is rows */
    215 void
    216 ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
    217                            CoinWorkDouble diagonalScaleFactor)
     213void ClpCholeskyBase::solveKKT(CoinWorkDouble *region1, CoinWorkDouble *region2, const CoinWorkDouble *diagonal,
     214  CoinWorkDouble diagonalScaleFactor)
    218215{
    219      if (!doKKT_) {
    220           int iColumn;
    221           int numberColumns = model_->numberColumns();
    222           int numberTotal = numberRows_ + numberColumns;
    223           CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal];
    224           for (iColumn = 0; iColumn < numberTotal; iColumn++) {
    225                region1[iColumn] *= diagonal[iColumn];
    226                region1Save[iColumn] = region1[iColumn];
    227           }
    228           multiplyAdd(region1 + numberColumns, numberRows_, -1.0, region2, 1.0);
    229           model_->clpMatrix()->times(1.0, region1, region2);
    230           CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
    231           CoinWorkDouble scale = 1.0;
    232           CoinWorkDouble unscale = 1.0;
    233           if (maximumRHS > 1.0e-30) {
    234                if (maximumRHS <= 0.5) {
    235                     CoinWorkDouble factor = 2.0;
    236                     while (maximumRHS <= 0.5) {
    237                          maximumRHS *= factor;
    238                          scale *= factor;
    239                     } /* endwhile */
    240                } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
    241                     CoinWorkDouble factor = 0.5;
    242                     while (maximumRHS >= 2.0) {
    243                          maximumRHS *= factor;
    244                          scale *= factor;
    245                     } /* endwhile */
    246                }
    247                unscale = diagonalScaleFactor / scale;
    248           } else {
    249                //effectively zero
    250                scale = 0.0;
    251                unscale = 0.0;
    252           }
    253           multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
    254           solve(region2);
    255           multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
    256           multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns, 0.0);
    257           CoinZeroN(region1, numberColumns);
    258           model_->clpMatrix()->transposeTimes(1.0, region2, region1);
    259           for (iColumn = 0; iColumn < numberTotal; iColumn++)
    260                region1[iColumn] = region1[iColumn] * diagonal[iColumn] - region1Save[iColumn];
    261           delete [] region1Save;
    262      } else {
    263           // KKT
    264           int numberRowsModel = model_->numberRows();
    265           int numberColumns = model_->numberColumns();
    266           int numberTotal = numberColumns + numberRowsModel;
    267           CoinWorkDouble * array = new CoinWorkDouble [numberRows_];
    268           CoinMemcpyN(region1, numberTotal, array);
    269           CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
    270           assert (numberRows_ >= numberRowsModel + numberTotal);
    271           solve(array);
    272           int iRow;
    273           for (iRow = 0; iRow < numberTotal; iRow++) {
    274                if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
    275                 COIN_DETAIL_PRINT(printf("row region1 %d dropped %g\n", iRow, array[iRow]));
    276                }
    277           }
    278           for (; iRow < numberRows_; iRow++) {
    279                if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
    280                 COIN_DETAIL_PRINT(printf("row region2 %d dropped %g\n", iRow, array[iRow]));
    281                }
    282           }
    283           CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
    284           CoinMemcpyN(array, numberTotal, region1);
    285           delete [] array;
    286      }
     216  if (!doKKT_) {
     217    int iColumn;
     218    int numberColumns = model_->numberColumns();
     219    int numberTotal = numberRows_ + numberColumns;
     220    CoinWorkDouble *region1Save = new CoinWorkDouble[numberTotal];
     221    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
     222      region1[iColumn] *= diagonal[iColumn];
     223      region1Save[iColumn] = region1[iColumn];
     224    }
     225    multiplyAdd(region1 + numberColumns, numberRows_, -1.0, region2, 1.0);
     226    model_->clpMatrix()->times(1.0, region1, region2);
     227    CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
     228    CoinWorkDouble scale = 1.0;
     229    CoinWorkDouble unscale = 1.0;
     230    if (maximumRHS > 1.0e-30) {
     231      if (maximumRHS <= 0.5) {
     232        CoinWorkDouble factor = 2.0;
     233        while (maximumRHS <= 0.5) {
     234          maximumRHS *= factor;
     235          scale *= factor;
     236        } /* endwhile */
     237      } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
     238        CoinWorkDouble factor = 0.5;
     239        while (maximumRHS >= 2.0) {
     240          maximumRHS *= factor;
     241          scale *= factor;
     242        } /* endwhile */
     243      }
     244      unscale = diagonalScaleFactor / scale;
     245    } else {
     246      //effectively zero
     247      scale = 0.0;
     248      unscale = 0.0;
     249    }
     250    multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
     251    solve(region2);
     252    multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
     253    multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns, 0.0);
     254    CoinZeroN(region1, numberColumns);
     255    model_->clpMatrix()->transposeTimes(1.0, region2, region1);
     256    for (iColumn = 0; iColumn < numberTotal; iColumn++)
     257      region1[iColumn] = region1[iColumn] * diagonal[iColumn] - region1Save[iColumn];
     258    delete[] region1Save;
     259  } else {
     260    // KKT
     261    int numberRowsModel = model_->numberRows();
     262    int numberColumns = model_->numberColumns();
     263    int numberTotal = numberColumns + numberRowsModel;
     264    CoinWorkDouble *array = new CoinWorkDouble[numberRows_];
     265    CoinMemcpyN(region1, numberTotal, array);
     266    CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
     267    assert(numberRows_ >= numberRowsModel + numberTotal);
     268    solve(array);
     269    int iRow;
     270    for (iRow = 0; iRow < numberTotal; iRow++) {
     271      if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
     272        COIN_DETAIL_PRINT(printf("row region1 %d dropped %g\n", iRow, array[iRow]));
     273      }
     274    }
     275    for (; iRow < numberRows_; iRow++) {
     276      if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
     277        COIN_DETAIL_PRINT(printf("row region2 %d dropped %g\n", iRow, array[iRow]));
     278      }
     279    }
     280    CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
     281    CoinMemcpyN(array, numberTotal, region1);
     282    delete[] array;
     283  }
    287284}
    288285//-------------------------------------------------------------------
    289286// Clone
    290287//-------------------------------------------------------------------
    291 ClpCholeskyBase * ClpCholeskyBase::clone() const
     288ClpCholeskyBase *ClpCholeskyBase::clone() const
    292289{
    293      return new ClpCholeskyBase(*this);
     290  return new ClpCholeskyBase(*this);
    294291}
    295292// Forms ADAT - returns nonzero if not enough memory
    296 int
    297 ClpCholeskyBase::preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT)
     293int ClpCholeskyBase::preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT)
    298294{
    299      delete rowCopy_;
    300      rowCopy_ = model_->clpMatrix()->reverseOrderedCopy();
    301      if (!doKKT) {
    302           numberRows_ = model_->numberRows();
    303           rowsDropped_ = new char [numberRows_];
    304           memset(rowsDropped_, 0, numberRows_);
    305           numberRowsDropped_ = 0;
    306           // Space for starts
    307           choleskyStart_ = new int[numberRows_+1];
    308           const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    309           const int * columnLength = model_->clpMatrix()->getVectorLengths();
    310           const int * row = model_->clpMatrix()->getIndices();
    311           const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    312           const int * rowLength = rowCopy_->getVectorLengths();
    313           const int * column = rowCopy_->getIndices();
    314           // We need two arrays for counts
    315           int * which = new int [numberRows_];
    316           int * used = new int[numberRows_+1];
    317           CoinZeroN(used, numberRows_);
    318           int iRow;
    319           sizeFactor_ = 0;
    320           int numberColumns = model_->numberColumns();
    321           int numberDense = 0;
    322           //denseThreshold_=3;
    323           if (denseThreshold_ > 0) {
    324                delete [] whichDense_;
    325                delete [] denseColumn_;
    326                delete dense_;
    327                whichDense_ = new char[numberColumns];
    328                int iColumn;
    329                used[numberRows_] = 0;
    330                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    331                     int length = columnLength[iColumn];
    332                     used[length] += 1;
    333                }
    334                int nLong = 0;
    335                int stop = CoinMax(denseThreshold_ / 2, 100);
    336                for (iRow = numberRows_; iRow >= stop; iRow--) {
    337                     if (used[iRow])
    338                       COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
    339                     nLong += used[iRow];
    340                     if (nLong > 50 || nLong > (numberColumns >> 2))
    341                          break;
    342                }
    343                CoinZeroN(used, numberRows_);
    344                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    345                     if (columnLength[iColumn] < denseThreshold_) {
    346                          whichDense_[iColumn] = 0;
    347                     } else {
    348                          whichDense_[iColumn] = 1;
    349                          numberDense++;
    350                     }
    351                }
    352                if (!numberDense || numberDense > 100) {
    353                     // free
    354                     delete [] whichDense_;
    355                     whichDense_ = NULL;
    356                     denseColumn_ = NULL;
    357                     dense_ = NULL;
    358                } else {
    359                     // space for dense columns
    360                     denseColumn_ = new longDouble [numberDense*numberRows_];
    361                     // dense cholesky
    362                     dense_ = new ClpCholeskyDense();
    363                     dense_->reserveSpace(NULL, numberDense);
    364                     COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
    365                }
    366           }
    367           int offset = includeDiagonal ? 0 : 1;
    368           if (lowerTriangular)
    369                offset = -offset;
    370           for (iRow = 0; iRow < numberRows_; iRow++) {
    371                int number = 0;
    372                // make sure diagonal exists if includeDiagonal
    373                if (!offset) {
    374                     which[0] = iRow;
    375                     used[iRow] = 1;
    376                     number = 1;
    377                }
    378                CoinBigIndex startRow = rowStart[iRow];
    379                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    380                if (lowerTriangular) {
    381                     for (CoinBigIndex k = startRow; k < endRow; k++) {
    382                          int iColumn = column[k];
    383                          if (!whichDense_ || !whichDense_[iColumn]) {
    384                               CoinBigIndex start = columnStart[iColumn];
    385                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    386                               for (CoinBigIndex j = start; j < end; j++) {
    387                                    int jRow = row[j];
    388                                    if (jRow <= iRow + offset) {
    389                                         if (!used[jRow]) {
    390                                              used[jRow] = 1;
    391                                              which[number++] = jRow;
    392                                         }
    393                                    }
    394                               }
    395                          }
    396                     }
    397                } else {
    398                     for (CoinBigIndex k = startRow; k < endRow; k++) {
    399                          int iColumn = column[k];
    400                          if (!whichDense_ || !whichDense_[iColumn]) {
    401                               CoinBigIndex start = columnStart[iColumn];
    402                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    403                               for (CoinBigIndex j = start; j < end; j++) {
    404                                    int jRow = row[j];
    405                                    if (jRow >= iRow + offset) {
    406                                         if (!used[jRow]) {
    407                                              used[jRow] = 1;
    408                                              which[number++] = jRow;
    409                                         }
    410                                    }
    411                               }
    412                          }
    413                     }
    414                }
    415                sizeFactor_ += number;
    416                int j;
    417                for (j = 0; j < number; j++)
    418                     used[which[j]] = 0;
    419           }
    420           delete [] which;
    421           // Now we have size - create arrays and fill in
    422           try {
    423                choleskyRow_ = new int [sizeFactor_];
    424           } catch (...) {
    425                // no memory
    426                delete [] choleskyStart_;
    427                choleskyStart_ = NULL;
    428                return -1;
    429           }
    430           sizeFactor_ = 0;
    431           which = choleskyRow_;
    432           for (iRow = 0; iRow < numberRows_; iRow++) {
    433                int number = 0;
    434                // make sure diagonal exists if includeDiagonal
    435                if (!offset) {
    436                     which[0] = iRow;
    437                     used[iRow] = 1;
    438                     number = 1;
    439                }
    440                choleskyStart_[iRow] = sizeFactor_;
    441                CoinBigIndex startRow = rowStart[iRow];
    442                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    443                if (lowerTriangular) {
    444                     for (CoinBigIndex k = startRow; k < endRow; k++) {
    445                          int iColumn = column[k];
    446                          if (!whichDense_ || !whichDense_[iColumn]) {
    447                               CoinBigIndex start = columnStart[iColumn];
    448                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    449                               for (CoinBigIndex j = start; j < end; j++) {
    450                                    int jRow = row[j];
    451                                    if (jRow <= iRow + offset) {
    452                                         if (!used[jRow]) {
    453                                              used[jRow] = 1;
    454                                              which[number++] = jRow;
    455                                         }
    456                                    }
    457                               }
    458                          }
    459                     }
    460                } else {
    461                     for (CoinBigIndex k = startRow; k < endRow; k++) {
    462                          int iColumn = column[k];
    463                          if (!whichDense_ || !whichDense_[iColumn]) {
    464                               CoinBigIndex start = columnStart[iColumn];
    465                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    466                               for (CoinBigIndex j = start; j < end; j++) {
    467                                    int jRow = row[j];
    468                                    if (jRow >= iRow + offset) {
    469                                         if (!used[jRow]) {
    470                                              used[jRow] = 1;
    471                                              which[number++] = jRow;
    472                                         }
    473                                    }
    474                               }
    475                          }
    476                     }
    477                }
    478                sizeFactor_ += number;
    479                int j;
    480                for (j = 0; j < number; j++)
    481                     used[which[j]] = 0;
    482                // Sort
    483                std::sort(which, which + number);
    484                // move which on
    485                which += number;
    486           }
    487           choleskyStart_[numberRows_] = sizeFactor_;
    488           delete [] used;
    489           return 0;
    490      } else {
    491           int numberRowsModel = model_->numberRows();
    492           int numberColumns = model_->numberColumns();
    493           int numberTotal = numberColumns + numberRowsModel;
    494           numberRows_ = 2 * numberRowsModel + numberColumns;
    495           rowsDropped_ = new char [numberRows_];
    496           memset(rowsDropped_, 0, numberRows_);
    497           numberRowsDropped_ = 0;
    498           CoinPackedMatrix * quadratic = NULL;
    499           ClpQuadraticObjective * quadraticObj =
    500                (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
    501           if (quadraticObj)
    502                quadratic = quadraticObj->quadraticObjective();
    503           CoinBigIndex numberElements = model_->clpMatrix()->getNumElements();
    504           numberElements = numberElements + 2 * numberRowsModel + numberTotal;
    505           if (quadratic)
    506                numberElements += quadratic->getNumElements();
    507           // Space for starts
    508           choleskyStart_ = new int[numberRows_+1];
    509           const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    510           const int * columnLength = model_->clpMatrix()->getVectorLengths();
    511           const int * row = model_->clpMatrix()->getIndices();
    512           //const double * element = model_->clpMatrix()->getElements();
    513           // Now we have size - create arrays and fill in
    514           try {
    515                choleskyRow_ = new int [numberElements];
    516           } catch (...) {
    517                // no memory
    518                delete [] choleskyStart_;
    519                choleskyStart_ = NULL;
    520                return -1;
    521           }
    522           int iRow, iColumn;
     295  delete rowCopy_;
     296  rowCopy_ = model_->clpMatrix()->reverseOrderedCopy();
     297  if (!doKKT) {
     298    numberRows_ = model_->numberRows();
     299    rowsDropped_ = new char[numberRows_];
     300    memset(rowsDropped_, 0, numberRows_);
     301    numberRowsDropped_ = 0;
     302    // Space for starts
     303    choleskyStart_ = new int[numberRows_ + 1];
     304    const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     305    const int *columnLength = model_->clpMatrix()->getVectorLengths();
     306    const int *row = model_->clpMatrix()->getIndices();
     307    const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     308    const int *rowLength = rowCopy_->getVectorLengths();
     309    const int *column = rowCopy_->getIndices();
     310    // We need two arrays for counts
     311    int *which = new int[numberRows_];
     312    int *used = new int[numberRows_ + 1];
     313    CoinZeroN(used, numberRows_);
     314    int iRow;
     315    sizeFactor_ = 0;
     316    int numberColumns = model_->numberColumns();
     317    int numberDense = 0;
     318    //denseThreshold_=3;
     319    if (denseThreshold_ > 0) {
     320      delete[] whichDense_;
     321      delete[] denseColumn_;
     322      delete dense_;
     323      whichDense_ = new char[numberColumns];
     324      int iColumn;
     325      used[numberRows_] = 0;
     326      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     327        int length = columnLength[iColumn];
     328        used[length] += 1;
     329      }
     330      int nLong = 0;
     331      int stop = CoinMax(denseThreshold_ / 2, 100);
     332      for (iRow = numberRows_; iRow >= stop; iRow--) {
     333        if (used[iRow])
     334          COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
     335        nLong += used[iRow];
     336        if (nLong > 50 || nLong > (numberColumns >> 2))
     337          break;
     338      }
     339      CoinZeroN(used, numberRows_);
     340      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     341        if (columnLength[iColumn] < denseThreshold_) {
     342          whichDense_[iColumn] = 0;
     343        } else {
     344          whichDense_[iColumn] = 1;
     345          numberDense++;
     346        }
     347      }
     348      if (!numberDense || numberDense > 100) {
     349        // free
     350        delete[] whichDense_;
     351        whichDense_ = NULL;
     352        denseColumn_ = NULL;
     353        dense_ = NULL;
     354      } else {
     355        // space for dense columns
     356        denseColumn_ = new longDouble[numberDense * numberRows_];
     357        // dense cholesky
     358        dense_ = new ClpCholeskyDense();
     359        dense_->reserveSpace(NULL, numberDense);
     360        COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
     361      }
     362    }
     363    int offset = includeDiagonal ? 0 : 1;
     364    if (lowerTriangular)
     365      offset = -offset;
     366    for (iRow = 0; iRow < numberRows_; iRow++) {
     367      int number = 0;
     368      // make sure diagonal exists if includeDiagonal
     369      if (!offset) {
     370        which[0] = iRow;
     371        used[iRow] = 1;
     372        number = 1;
     373      }
     374      CoinBigIndex startRow = rowStart[iRow];
     375      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     376      if (lowerTriangular) {
     377        for (CoinBigIndex k = startRow; k < endRow; k++) {
     378          int iColumn = column[k];
     379          if (!whichDense_ || !whichDense_[iColumn]) {
     380            CoinBigIndex start = columnStart[iColumn];
     381            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     382            for (CoinBigIndex j = start; j < end; j++) {
     383              int jRow = row[j];
     384              if (jRow <= iRow + offset) {
     385                if (!used[jRow]) {
     386                  used[jRow] = 1;
     387                  which[number++] = jRow;
     388                }
     389              }
     390            }
     391          }
     392        }
     393      } else {
     394        for (CoinBigIndex k = startRow; k < endRow; k++) {
     395          int iColumn = column[k];
     396          if (!whichDense_ || !whichDense_[iColumn]) {
     397            CoinBigIndex start = columnStart[iColumn];
     398            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     399            for (CoinBigIndex j = start; j < end; j++) {
     400              int jRow = row[j];
     401              if (jRow >= iRow + offset) {
     402                if (!used[jRow]) {
     403                  used[jRow] = 1;
     404                  which[number++] = jRow;
     405                }
     406              }
     407            }
     408          }
     409        }
     410      }
     411      sizeFactor_ += number;
     412      int j;
     413      for (j = 0; j < number; j++)
     414        used[which[j]] = 0;
     415    }
     416    delete[] which;
     417    // Now we have size - create arrays and fill in
     418    try {
     419      choleskyRow_ = new int[sizeFactor_];
     420    } catch (...) {
     421      // no memory
     422      delete[] choleskyStart_;
     423      choleskyStart_ = NULL;
     424      return -1;
     425    }
     426    sizeFactor_ = 0;
     427    which = choleskyRow_;
     428    for (iRow = 0; iRow < numberRows_; iRow++) {
     429      int number = 0;
     430      // make sure diagonal exists if includeDiagonal
     431      if (!offset) {
     432        which[0] = iRow;
     433        used[iRow] = 1;
     434        number = 1;
     435      }
     436      choleskyStart_[iRow] = sizeFactor_;
     437      CoinBigIndex startRow = rowStart[iRow];
     438      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     439      if (lowerTriangular) {
     440        for (CoinBigIndex k = startRow; k < endRow; k++) {
     441          int iColumn = column[k];
     442          if (!whichDense_ || !whichDense_[iColumn]) {
     443            CoinBigIndex start = columnStart[iColumn];
     444            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     445            for (CoinBigIndex j = start; j < end; j++) {
     446              int jRow = row[j];
     447              if (jRow <= iRow + offset) {
     448                if (!used[jRow]) {
     449                  used[jRow] = 1;
     450                  which[number++] = jRow;
     451                }
     452              }
     453            }
     454          }
     455        }
     456      } else {
     457        for (CoinBigIndex k = startRow; k < endRow; k++) {
     458          int iColumn = column[k];
     459          if (!whichDense_ || !whichDense_[iColumn]) {
     460            CoinBigIndex start = columnStart[iColumn];
     461            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     462            for (CoinBigIndex j = start; j < end; j++) {
     463              int jRow = row[j];
     464              if (jRow >= iRow + offset) {
     465                if (!used[jRow]) {
     466                  used[jRow] = 1;
     467                  which[number++] = jRow;
     468                }
     469              }
     470            }
     471          }
     472        }
     473      }
     474      sizeFactor_ += number;
     475      int j;
     476      for (j = 0; j < number; j++)
     477        used[which[j]] = 0;
     478      // Sort
     479      std::sort(which, which + number);
     480      // move which on
     481      which += number;
     482    }
     483    choleskyStart_[numberRows_] = sizeFactor_;
     484    delete[] used;
     485    return 0;
     486  } else {
     487    int numberRowsModel = model_->numberRows();
     488    int numberColumns = model_->numberColumns();
     489    int numberTotal = numberColumns + numberRowsModel;
     490    numberRows_ = 2 * numberRowsModel + numberColumns;
     491    rowsDropped_ = new char[numberRows_];
     492    memset(rowsDropped_, 0, numberRows_);
     493    numberRowsDropped_ = 0;
     494    CoinPackedMatrix *quadratic = NULL;
     495    ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
     496    if (quadraticObj)
     497      quadratic = quadraticObj->quadraticObjective();
     498    CoinBigIndex numberElements = model_->clpMatrix()->getNumElements();
     499    numberElements = numberElements + 2 * numberRowsModel + numberTotal;
     500    if (quadratic)
     501      numberElements += quadratic->getNumElements();
     502    // Space for starts
     503    choleskyStart_ = new int[numberRows_ + 1];
     504    const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     505    const int *columnLength = model_->clpMatrix()->getVectorLengths();
     506    const int *row = model_->clpMatrix()->getIndices();
     507    //const double * element = model_->clpMatrix()->getElements();
     508    // Now we have size - create arrays and fill in
     509    try {
     510      choleskyRow_ = new int[numberElements];
     511    } catch (...) {
     512      // no memory
     513      delete[] choleskyStart_;
     514      choleskyStart_ = NULL;
     515      return -1;
     516    }
     517    int iRow, iColumn;
    523518
    524           sizeFactor_ = 0;
    525           // matrix
    526           if (lowerTriangular) {
    527                if (!quadratic) {
    528                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    529                          choleskyStart_[iColumn] = sizeFactor_;
    530                          choleskyRow_[sizeFactor_++] = iColumn;
    531                          CoinBigIndex start = columnStart[iColumn];
    532                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    533                          if (!includeDiagonal)
    534                               start++;
    535                          for (CoinBigIndex j = start; j < end; j++) {
    536                               choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
    537                          }
    538                     }
    539                } else {
    540                     // Quadratic
    541                     const int * columnQuadratic = quadratic->getIndices();
    542                     const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    543                     const int * columnQuadraticLength = quadratic->getVectorLengths();
    544                     //const double * quadraticElement = quadratic->getElements();
    545                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    546                          choleskyStart_[iColumn] = sizeFactor_;
    547                          if (includeDiagonal)
    548                               choleskyRow_[sizeFactor_++] = iColumn;
    549                          CoinBigIndex j;
    550                          for ( j = columnQuadraticStart[iColumn];
    551                                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    552                               int jColumn = columnQuadratic[j];
    553                               if (jColumn > iColumn)
    554                                    choleskyRow_[sizeFactor_++] = jColumn;
    555                          }
    556                          CoinBigIndex start = columnStart[iColumn];
    557                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    558                          for ( j = start; j < end; j++) {
    559                               choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
    560                          }
    561                     }
    562                }
    563                // slacks
    564                for (; iColumn < numberTotal; iColumn++) {
    565                     choleskyStart_[iColumn] = sizeFactor_;
    566                     if (includeDiagonal)
    567                          choleskyRow_[sizeFactor_++] = iColumn;
    568                     choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
    569                }
    570                // Transpose - nonzero diagonal (may regularize)
    571                for (iRow = 0; iRow < numberRowsModel; iRow++) {
    572                     choleskyStart_[iRow+numberTotal] = sizeFactor_;
    573                     // diagonal
    574                     if (includeDiagonal)
    575                          choleskyRow_[sizeFactor_++] = iRow + numberTotal;
    576                }
    577                choleskyStart_[numberRows_] = sizeFactor_;
    578           } else {
    579                // transpose
    580                ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
    581                const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    582                const int * rowLength = rowCopy->getVectorLengths();
    583                const int * column = rowCopy->getIndices();
    584                if (!quadratic) {
    585                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    586                          choleskyStart_[iColumn] = sizeFactor_;
    587                          if (includeDiagonal)
    588                               choleskyRow_[sizeFactor_++] = iColumn;
    589                     }
    590                } else {
    591                     // Quadratic
    592                     // transpose
    593                     CoinPackedMatrix quadraticT;
    594                     quadraticT.reverseOrderedCopyOf(*quadratic);
    595                     const int * columnQuadratic = quadraticT.getIndices();
    596                     const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
    597                     const int * columnQuadraticLength = quadraticT.getVectorLengths();
    598                     //const double * quadraticElement = quadraticT.getElements();
    599                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    600                          choleskyStart_[iColumn] = sizeFactor_;
    601                          for (CoinBigIndex j = columnQuadraticStart[iColumn];
    602                                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    603                               int jColumn = columnQuadratic[j];
    604                               if (jColumn < iColumn)
    605                                    choleskyRow_[sizeFactor_++] = jColumn;
    606                          }
    607                          if (includeDiagonal)
    608                               choleskyRow_[sizeFactor_++] = iColumn;
    609                     }
    610                }
    611                int iRow;
    612                // slacks
    613                for (iRow = 0; iRow < numberRowsModel; iRow++) {
    614                     choleskyStart_[iRow+numberColumns] = sizeFactor_;
    615                     if (includeDiagonal)
    616                          choleskyRow_[sizeFactor_++] = iRow + numberColumns;
    617                }
    618                for (iRow = 0; iRow < numberRowsModel; iRow++) {
    619                     choleskyStart_[iRow+numberTotal] = sizeFactor_;
    620                     CoinBigIndex start = rowStart[iRow];
    621                     CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
    622                     for (CoinBigIndex j = start; j < end; j++) {
    623                          choleskyRow_[sizeFactor_++] = column[j];
    624                     }
    625                     // slack
    626                     choleskyRow_[sizeFactor_++] = numberColumns + iRow;
    627                     if (includeDiagonal)
    628                          choleskyRow_[sizeFactor_++] = iRow + numberTotal;
    629                }
    630                choleskyStart_[numberRows_] = sizeFactor_;
    631           }
    632      }
    633      return 0;
     519    sizeFactor_ = 0;
     520    // matrix
     521    if (lowerTriangular) {
     522      if (!quadratic) {
     523        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     524          choleskyStart_[iColumn] = sizeFactor_;
     525          choleskyRow_[sizeFactor_++] = iColumn;
     526          CoinBigIndex start = columnStart[iColumn];
     527          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     528          if (!includeDiagonal)
     529            start++;
     530          for (CoinBigIndex j = start; j < end; j++) {
     531            choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
     532          }
     533        }
     534      } else {
     535        // Quadratic
     536        const int *columnQuadratic = quadratic->getIndices();
     537        const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
     538        const int *columnQuadraticLength = quadratic->getVectorLengths();
     539        //const double * quadraticElement = quadratic->getElements();
     540        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     541          choleskyStart_[iColumn] = sizeFactor_;
     542          if (includeDiagonal)
     543            choleskyRow_[sizeFactor_++] = iColumn;
     544          CoinBigIndex j;
     545          for (j = columnQuadraticStart[iColumn];
     546               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     547            int jColumn = columnQuadratic[j];
     548            if (jColumn > iColumn)
     549              choleskyRow_[sizeFactor_++] = jColumn;
     550          }
     551          CoinBigIndex start = columnStart[iColumn];
     552          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     553          for (j = start; j < end; j++) {
     554            choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
     555          }
     556        }
     557      }
     558      // slacks
     559      for (; iColumn < numberTotal; iColumn++) {
     560        choleskyStart_[iColumn] = sizeFactor_;
     561        if (includeDiagonal)
     562          choleskyRow_[sizeFactor_++] = iColumn;
     563        choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
     564      }
     565      // Transpose - nonzero diagonal (may regularize)
     566      for (iRow = 0; iRow < numberRowsModel; iRow++) {
     567        choleskyStart_[iRow + numberTotal] = sizeFactor_;
     568        // diagonal
     569        if (includeDiagonal)
     570          choleskyRow_[sizeFactor_++] = iRow + numberTotal;
     571      }
     572      choleskyStart_[numberRows_] = sizeFactor_;
     573    } else {
     574      // transpose
     575      ClpMatrixBase *rowCopy = model_->clpMatrix()->reverseOrderedCopy();
     576      const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
     577      const int *rowLength = rowCopy->getVectorLengths();
     578      const int *column = rowCopy->getIndices();
     579      if (!quadratic) {
     580        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     581          choleskyStart_[iColumn] = sizeFactor_;
     582          if (includeDiagonal)
     583            choleskyRow_[sizeFactor_++] = iColumn;
     584        }
     585      } else {
     586        // Quadratic
     587        // transpose
     588        CoinPackedMatrix quadraticT;
     589        quadraticT.reverseOrderedCopyOf(*quadratic);
     590        const int *columnQuadratic = quadraticT.getIndices();
     591        const CoinBigIndex *columnQuadraticStart = quadraticT.getVectorStarts();
     592        const int *columnQuadraticLength = quadraticT.getVectorLengths();
     593        //const double * quadraticElement = quadraticT.getElements();
     594        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     595          choleskyStart_[iColumn] = sizeFactor_;
     596          for (CoinBigIndex j = columnQuadraticStart[iColumn];
     597               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     598            int jColumn = columnQuadratic[j];
     599            if (jColumn < iColumn)
     600              choleskyRow_[sizeFactor_++] = jColumn;
     601          }
     602          if (includeDiagonal)
     603            choleskyRow_[sizeFactor_++] = iColumn;
     604        }
     605      }
     606      int iRow;
     607      // slacks
     608      for (iRow = 0; iRow < numberRowsModel; iRow++) {
     609        choleskyStart_[iRow + numberColumns] = sizeFactor_;
     610        if (includeDiagonal)
     611          choleskyRow_[sizeFactor_++] = iRow + numberColumns;
     612      }
     613      for (iRow = 0; iRow < numberRowsModel; iRow++) {
     614        choleskyStart_[iRow + numberTotal] = sizeFactor_;
     615        CoinBigIndex start = rowStart[iRow];
     616        CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
     617        for (CoinBigIndex j = start; j < end; j++) {
     618          choleskyRow_[sizeFactor_++] = column[j];
     619        }
     620        // slack
     621        choleskyRow_[sizeFactor_++] = numberColumns + iRow;
     622        if (includeDiagonal)
     623          choleskyRow_[sizeFactor_++] = iRow + numberTotal;
     624      }
     625      choleskyStart_[numberRows_] = sizeFactor_;
     626    }
     627  }
     628  return 0;
    634629}
    635630/* Orders rows and saves pointer to matrix.and model */
    636 int
    637 ClpCholeskyBase::order(ClpInterior * model)
     631int ClpCholeskyBase::order(ClpInterior *model)
    638632{
    639      model_ = model;
     633  model_ = model;
    640634#define BASE_ORDER 2
    641 #if BASE_ORDER>0
    642      if (model_->numberRows() > 6 ) {
    643           if (preOrder(doKKT_, true, doKKT_))
    644                return -1;
    645           //rowsDropped_ = new char [numberRows_];
    646           numberRowsDropped_ = 0;
    647           memset(rowsDropped_, 0, numberRows_);
    648           //rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
    649           // approximate minimum degree
    650           return orderAMD();
    651      }
    652 #endif
    653      int numberRowsModel = model_->numberRows();
    654      int numberColumns = model_->numberColumns();
    655      int numberTotal = numberColumns + numberRowsModel;
    656      CoinPackedMatrix * quadratic = NULL;
    657      ClpQuadraticObjective * quadraticObj =
    658           (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
    659      if (quadraticObj)
    660           quadratic = quadraticObj->quadraticObjective();
    661      if (!doKKT_) {
    662           numberRows_ = model->numberRows();
    663      } else {
    664           numberRows_ = 2 * numberRowsModel + numberColumns;
    665      }
    666      rowsDropped_ = new char [numberRows_];
    667      numberRowsDropped_ = 0;
    668      memset(rowsDropped_, 0, numberRows_);
    669      rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
    670      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    671      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    672      const int * row = model_->clpMatrix()->getIndices();
    673      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    674      const int * rowLength = rowCopy_->getVectorLengths();
    675      const int * column = rowCopy_->getIndices();
    676      // We need arrays for counts
    677      int * which = new int [numberRows_];
    678      int * used = new int[numberRows_+1];
    679      int * count = new int[numberRows_];
    680      CoinZeroN(count, numberRows_);
    681      CoinZeroN(used, numberRows_);
    682      int iRow;
    683      sizeFactor_ = 0;
    684      permute_ = new int[numberRows_];
    685      for (iRow = 0; iRow < numberRows_; iRow++)
    686           permute_[iRow] = iRow;
    687      if (!doKKT_) {
    688           int numberDense = 0;
    689           if (denseThreshold_ > 0) {
    690                delete [] whichDense_;
    691                delete [] denseColumn_;
    692                delete dense_;
    693                whichDense_ = new char[numberColumns];
    694                int iColumn;
    695                used[numberRows_] = 0;
    696                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    697                     int length = columnLength[iColumn];
    698                     used[length] += 1;
    699                }
    700                int nLong = 0;
    701                int stop = CoinMax(denseThreshold_ / 2, 100);
    702                for (iRow = numberRows_; iRow >= stop; iRow--) {
    703                     if (used[iRow])
    704                       COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
    705                     nLong += used[iRow];
    706                     if (nLong > 50 || nLong > (numberColumns >> 2))
    707                          break;
    708                }
    709                CoinZeroN(used, numberRows_);
    710                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    711                     if (columnLength[iColumn] < denseThreshold_) {
    712                          whichDense_[iColumn] = 0;
    713                     } else {
    714                          whichDense_[iColumn] = 1;
    715                          numberDense++;
    716                     }
    717                }
    718                if (!numberDense || numberDense > 100) {
    719                     // free
    720                     delete [] whichDense_;
    721                     whichDense_ = NULL;
    722                     denseColumn_ = NULL;
    723                     dense_ = NULL;
    724                } else {
    725                     // space for dense columns
    726                     denseColumn_ = new longDouble [numberDense*numberRows_];
    727                     // dense cholesky
    728                     dense_ = new ClpCholeskyDense();
    729                     dense_->reserveSpace(NULL, numberDense);
    730                     COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
    731                }
    732           }
    733           /*
     635#if BASE_ORDER > 0
     636  if (model_->numberRows() > 6) {
     637    if (preOrder(doKKT_, true, doKKT_))
     638      return -1;
     639    //rowsDropped_ = new char [numberRows_];
     640    numberRowsDropped_ = 0;
     641    memset(rowsDropped_, 0, numberRows_);
     642    //rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     643    // approximate minimum degree
     644    return orderAMD();
     645  }
     646#endif
     647  int numberRowsModel = model_->numberRows();
     648  int numberColumns = model_->numberColumns();
     649  int numberTotal = numberColumns + numberRowsModel;
     650  CoinPackedMatrix *quadratic = NULL;
     651  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
     652  if (quadraticObj)
     653    quadratic = quadraticObj->quadraticObjective();
     654  if (!doKKT_) {
     655    numberRows_ = model->numberRows();
     656  } else {
     657    numberRows_ = 2 * numberRowsModel + numberColumns;
     658  }
     659  rowsDropped_ = new char[numberRows_];
     660  numberRowsDropped_ = 0;
     661  memset(rowsDropped_, 0, numberRows_);
     662  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     663  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     664  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     665  const int *row = model_->clpMatrix()->getIndices();
     666  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     667  const int *rowLength = rowCopy_->getVectorLengths();
     668  const int *column = rowCopy_->getIndices();
     669  // We need arrays for counts
     670  int *which = new int[numberRows_];
     671  int *used = new int[numberRows_ + 1];
     672  int *count = new int[numberRows_];
     673  CoinZeroN(count, numberRows_);
     674  CoinZeroN(used, numberRows_);
     675  int iRow;
     676  sizeFactor_ = 0;
     677  permute_ = new int[numberRows_];
     678  for (iRow = 0; iRow < numberRows_; iRow++)
     679    permute_[iRow] = iRow;
     680  if (!doKKT_) {
     681    int numberDense = 0;
     682    if (denseThreshold_ > 0) {
     683      delete[] whichDense_;
     684      delete[] denseColumn_;
     685      delete dense_;
     686      whichDense_ = new char[numberColumns];
     687      int iColumn;
     688      used[numberRows_] = 0;
     689      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     690        int length = columnLength[iColumn];
     691        used[length] += 1;
     692      }
     693      int nLong = 0;
     694      int stop = CoinMax(denseThreshold_ / 2, 100);
     695      for (iRow = numberRows_; iRow >= stop; iRow--) {
     696        if (used[iRow])
     697          COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
     698        nLong += used[iRow];
     699        if (nLong > 50 || nLong > (numberColumns >> 2))
     700          break;
     701      }
     702      CoinZeroN(used, numberRows_);
     703      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     704        if (columnLength[iColumn] < denseThreshold_) {
     705          whichDense_[iColumn] = 0;
     706        } else {
     707          whichDense_[iColumn] = 1;
     708          numberDense++;
     709        }
     710      }
     711      if (!numberDense || numberDense > 100) {
     712        // free
     713        delete[] whichDense_;
     714        whichDense_ = NULL;
     715        denseColumn_ = NULL;
     716        dense_ = NULL;
     717      } else {
     718        // space for dense columns
     719        denseColumn_ = new longDouble[numberDense * numberRows_];
     720        // dense cholesky
     721        dense_ = new ClpCholeskyDense();
     722        dense_->reserveSpace(NULL, numberDense);
     723        COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
     724      }
     725    }
     726    /*
    734727             Get row counts and size
    735728          */
    736           for (iRow = 0; iRow < numberRows_; iRow++) {
    737                int number = 1;
    738                // make sure diagonal exists
    739                which[0] = iRow;
    740                used[iRow] = 1;
    741                CoinBigIndex startRow = rowStart[iRow];
    742                CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
    743                for (CoinBigIndex k = startRow; k < endRow; k++) {
    744                     int iColumn = column[k];
    745                     if (!whichDense_ || !whichDense_[iColumn]) {
    746                          CoinBigIndex start = columnStart[iColumn];
    747                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    748                          for (CoinBigIndex j = start; j < end; j++) {
    749                               int jRow = row[j];
    750                               if (jRow < iRow) {
    751                                    if (!used[jRow]) {
    752                                         used[jRow] = 1;
    753                                         which[number++] = jRow;
    754                                         count[jRow]++;
    755                                    }
    756                               }
    757                          }
    758                     }
    759                }
    760                sizeFactor_ += number;
    761                count[iRow] += number;
    762                int j;
    763                for (j = 0; j < number; j++)
    764                     used[which[j]] = 0;
    765           }
    766           CoinSort_2(count, count + numberRows_, permute_);
    767      } else {
    768           // KKT
    769           CoinBigIndex numberElements = model_->clpMatrix()->getNumElements();
    770           numberElements = numberElements + 2 * numberRowsModel + numberTotal;
    771           if (quadratic)
    772                numberElements += quadratic->getNumElements();
    773           // off diagonal
    774           numberElements -= numberRows_;
    775           sizeFactor_ = static_cast<int>(numberElements);
    776           // If we sort we need to redo symbolic etc
    777      }
    778      delete [] which;
    779      delete [] used;
    780      delete [] count;
    781      permuteInverse_ = new int [numberRows_];
    782      for (iRow = 0; iRow < numberRows_; iRow++) {
    783           //permute_[iRow]=iRow; // force no permute
    784           //permute_[iRow]=numberRows_-1-iRow; // force odd permute
    785           //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
    786           permuteInverse_[permute_[iRow]] = iRow;
    787      }
    788      return 0;
     729    for (iRow = 0; iRow < numberRows_; iRow++) {
     730      int number = 1;
     731      // make sure diagonal exists
     732      which[0] = iRow;
     733      used[iRow] = 1;
     734      CoinBigIndex startRow = rowStart[iRow];
     735      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
     736      for (CoinBigIndex k = startRow; k < endRow; k++) {
     737        int iColumn = column[k];
     738        if (!whichDense_ || !whichDense_[iColumn]) {
     739          CoinBigIndex start = columnStart[iColumn];
     740          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     741          for (CoinBigIndex j = start; j < end; j++) {
     742            int jRow = row[j];
     743            if (jRow < iRow) {
     744              if (!used[jRow]) {
     745                used[jRow] = 1;
     746                which[number++] = jRow;
     747                count[jRow]++;
     748              }
     749            }
     750          }
     751        }
     752      }
     753      sizeFactor_ += number;
     754      count[iRow] += number;
     755      int j;
     756      for (j = 0; j < number; j++)
     757        used[which[j]] = 0;
     758    }
     759    CoinSort_2(count, count + numberRows_, permute_);
     760  } else {
     761    // KKT
     762    CoinBigIndex numberElements = model_->clpMatrix()->getNumElements();
     763    numberElements = numberElements + 2 * numberRowsModel + numberTotal;
     764    if (quadratic)
     765      numberElements += quadratic->getNumElements();
     766    // off diagonal
     767    numberElements -= numberRows_;
     768    sizeFactor_ = static_cast< int >(numberElements);
     769    // If we sort we need to redo symbolic etc
     770  }
     771  delete[] which;
     772  delete[] used;
     773  delete[] count;
     774  permuteInverse_ = new int[numberRows_];
     775  for (iRow = 0; iRow < numberRows_; iRow++) {
     776    //permute_[iRow]=iRow; // force no permute
     777    //permute_[iRow]=numberRows_-1-iRow; // force odd permute
     778    //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
     779    permuteInverse_[permute_[iRow]] = iRow;
     780  }
     781  return 0;
    789782}
    790 #if BASE_ORDER==1
     783#if BASE_ORDER == 1
    791784/* Orders rows and saves pointer to matrix.and model */
    792 int
    793 ClpCholeskyBase::orderAMD()
     785int ClpCholeskyBase::orderAMD()
    794786{
    795      permuteInverse_ = new int [numberRows_];
    796      permute_ = new int[numberRows_];
    797      // Do ordering
    798      int returnCode = 0;
    799      // get more space and full matrix
    800      int space = 6 * sizeFactor_ + 100000;
    801      int * temp = new int [space];
    802      int * which = new int[2*numberRows_];
    803      CoinBigIndex * tempStart = new CoinBigIndex [numberRows_+1];
    804      memset(which, 0, numberRows_ * sizeof(int));
    805      for (int iRow = 0; iRow < numberRows_; iRow++) {
    806           which[iRow] += choleskyStart_[iRow+1] - choleskyStart_[iRow] - 1;
    807           for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
    808                int jRow = choleskyRow_[j];
    809                which[jRow]++;
    810           }
    811      }
    812      CoinBigIndex sizeFactor = 0;
    813      for (int iRow = 0; iRow < numberRows_; iRow++) {
    814           int length = which[iRow];
    815           permute_[iRow] = length;
    816           tempStart[iRow] = sizeFactor;
    817           which[iRow] = sizeFactor;
    818           sizeFactor += length;
    819      }
    820      for (int iRow = 0; iRow < numberRows_; iRow++) {
    821           assert (choleskyRow_[choleskyStart_[iRow]] == iRow);
    822           for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
    823                int jRow = choleskyRow_[j];
    824                int put = which[iRow];
    825                temp[put] = jRow;
    826                which[iRow]++;
    827                put = which[jRow];
    828                temp[put] = iRow;
    829                which[jRow]++;
    830           }
    831      }
    832      for (int iRow = 1; iRow < numberRows_; iRow++)
    833           assert (which[iRow-1] == tempStart[iRow]);
    834      CoinBigIndex lastSpace = sizeFactor;
    835      delete [] choleskyRow_;
    836      choleskyRow_ = temp;
    837      delete [] choleskyStart_;
    838      choleskyStart_ = tempStart;
    839      // linked lists of sizes and lengths
    840      int * first = new int [numberRows_];
    841      int * next = new int [numberRows_];
    842      int * previous = new int [numberRows_];
    843      char * mark = new char[numberRows_];
    844      memset(mark, 0, numberRows_);
    845      CoinBigIndex * sort = new CoinBigIndex [numberRows_];
    846      int * order = new int [numberRows_];
    847      for (int iRow = 0; iRow < numberRows_; iRow++) {
    848           first[iRow] = -1;
    849           next[iRow] = -1;
    850           previous[iRow] = -1;
    851           permuteInverse_[iRow] = -1;
    852      }
    853      int large = 1000 + 2 * numberRows_;
    854      for (int iRow = 0; iRow < numberRows_; iRow++) {
    855           int n = permute_[iRow];
    856           if (first[n] < 0) {
    857                first[n] = iRow;
    858                previous[iRow] = n + large;
    859                next[iRow] = n + 2 * large;
    860           } else {
    861                int k = first[n];
    862                assert (k < numberRows_);
    863                first[n] = iRow;
    864                previous[iRow] = n + large;
    865                assert (previous[k] == n + large);
    866                next[iRow] = k;
    867                previous[k] = iRow;
    868           }
    869      }
    870      int smallest = 0;
    871      int done = 0;
    872      int numberCompressions = 0;
    873      int numberExpansions = 0;
    874      while (smallest < numberRows_) {
    875           if (first[smallest] < 0 || first[smallest] > numberRows_) {
    876                smallest++;
    877                continue;
    878           }
    879           int iRow = first[smallest];
    880           if (false) {
    881                int kRow = -1;
    882                int ss = 999999;
    883                for (int jRow = numberRows_ - 1; jRow >= 0; jRow--) {
    884                     if (permuteInverse_[jRow] < 0) {
    885                          int length = permute_[jRow];
    886                          assert (length > 0);
    887                          for (CoinBigIndex j = choleskyStart_[jRow];
    888                                    j < choleskyStart_[jRow] + length; j++) {
    889                               int jjRow = choleskyRow_[j];
    890                               assert (permuteInverse_[jjRow] < 0);
    891                          }
    892                          if (length < ss) {
    893                               ss = length;
    894                               kRow = jRow;
    895                          }
    896                     }
    897                }
    898                assert (smallest == ss);
    899                printf("list chose %d - full chose %d - length %d\n",
    900                       iRow, kRow, ss);
    901           }
    902           int kNext = next[iRow];
    903           first[smallest] = kNext;
    904           if (kNext < numberRows_)
    905                previous[kNext] = previous[iRow];
    906           previous[iRow] = -1;
    907           next[iRow] = -1;
    908           permuteInverse_[iRow] = done;
    909           done++;
    910           // Now add edges
    911           CoinBigIndex start = choleskyStart_[iRow];
    912           CoinBigIndex end = choleskyStart_[iRow] + permute_[iRow];
    913           int nSave = 0;
    914           for (CoinBigIndex k = start; k < end; k++) {
    915                int kRow = choleskyRow_[k];
    916                assert (permuteInverse_[kRow] < 0);
    917                //if (permuteInverse_[kRow]<0)
    918                which[nSave++] = kRow;
    919           }
    920           for (int i = 0; i < nSave; i++) {
    921                int kRow = which[i];
    922                int length = permute_[kRow];
    923                CoinBigIndex start = choleskyStart_[kRow];
    924                CoinBigIndex end = choleskyStart_[kRow] + length;
    925                for (CoinBigIndex j = start; j < end; j++) {
    926                     int jRow = choleskyRow_[j];
    927                     mark[jRow] = 1;
    928                }
    929                mark[kRow] = 1;
    930                int n = nSave;
    931                for (int j = 0; j < nSave; j++) {
    932                     int jRow = which[j];
    933                     if (!mark[jRow]) {
    934                          which[n++] = jRow;
    935                     }
    936                }
    937                for (CoinBigIndex j = start; j < end; j++) {
    938                     int jRow = choleskyRow_[j];
    939                     mark[jRow] = 0;
    940                }
    941                mark[kRow] = 0;
    942                if (n > nSave) {
    943                     bool inPlace = (n - nSave == 1);
    944                     if (!inPlace) {
    945                          // extend
    946                          int length = n - nSave + end - start;
    947                          if (length + lastSpace > space) {
    948                               // need to compress
    949                               numberCompressions++;
    950                               int newN = 0;
    951                               for (int iRow = 0; iRow < numberRows_; iRow++) {
    952                                    CoinBigIndex start = choleskyStart_[iRow];
    953                                    if (permuteInverse_[iRow] < 0) {
    954                                         sort[newN] = start;
    955                                         order[newN++] = iRow;
    956                                    } else {
    957                                         choleskyStart_[iRow] = -1;
    958                                         permute_[iRow] = 0;
    959                                    }
    960                               }
    961                               CoinSort_2(sort, sort + newN, order);
    962                               sizeFactor = 0;
    963                               for (int k = 0; k < newN; k++) {
    964                                    int iRow = order[k];
    965                                    int length = permute_[iRow];
    966                                    CoinBigIndex start = choleskyStart_[iRow];
    967                                    choleskyStart_[iRow] = sizeFactor;
    968                                    for (int j = 0; j < length; j++)
    969                                         choleskyRow_[sizeFactor+j] = choleskyRow_[start+j];
    970                                    sizeFactor += length;
    971                               }
    972                               lastSpace = sizeFactor;
    973                               if ((sizeFactor + numberRows_ + 1000) * 4 > 3 * space) {
    974                                    // need to expand
    975                                    numberExpansions++;
    976                                    space = (3 * space) / 2;
    977                                    int * temp = new int [space];
    978                                    memcpy(temp, choleskyRow_, sizeFactor * sizeof(int));
    979                                    delete [] choleskyRow_;
    980                                    choleskyRow_ = temp;
    981                               }
    982                          }
    983                     }
    984                     // Now add
    985                     start = choleskyStart_[kRow];
    986                     end = choleskyStart_[kRow] + permute_[kRow];
    987                     if (!inPlace)
    988                          choleskyStart_[kRow] = lastSpace;
    989                     CoinBigIndex put = choleskyStart_[kRow];
    990                     for (CoinBigIndex j = start; j < end; j++) {
    991                          int jRow = choleskyRow_[j];
    992                          if (permuteInverse_[jRow] < 0)
    993                               choleskyRow_[put++] = jRow;
    994                     }
    995                     for (int j = nSave; j < n; j++) {
    996                          int jRow = which[j];
    997                          choleskyRow_[put++] = jRow;
    998                     }
    999                     if (!inPlace) {
    1000                          permute_[kRow] = put - lastSpace;
    1001                          lastSpace = put;
    1002                     } else {
    1003                          permute_[kRow] = put - choleskyStart_[kRow];
    1004                     }
    1005                } else {
    1006                     // pack down for new counts
    1007                     CoinBigIndex put = start;
    1008                     for (CoinBigIndex j = start; j < end; j++) {
    1009                          int jRow = choleskyRow_[j];
    1010                          if (permuteInverse_[jRow] < 0)
    1011                               choleskyRow_[put++] = jRow;
    1012                     }
    1013                     permute_[kRow] = put - start;
    1014                }
    1015                // take out
    1016                int iNext = next[kRow];
    1017                int iPrevious = previous[kRow];
    1018                if (iPrevious < numberRows_) {
    1019                     next[iPrevious] = iNext;
    1020                } else {
    1021                     assert (iPrevious == length + large);
    1022                     first[length] = iNext;
    1023                }
    1024                if (iNext < numberRows_) {
    1025                     previous[iNext] = iPrevious;
    1026                } else {
    1027                     assert (iNext == length + 2 * large);
    1028                }
    1029                // put in
    1030                length = permute_[kRow];
    1031                smallest = CoinMin(smallest, length);
    1032                if (first[length] < 0 || first[length] > numberRows_) {
    1033                     first[length] = kRow;
    1034                     previous[kRow] = length + large;
    1035                     next[kRow] = length + 2 * large;
    1036                } else {
    1037                     int k = first[length];
    1038                     assert (k < numberRows_);
    1039                     first[length] = kRow;
    1040                     previous[kRow] = length + large;
    1041                     assert (previous[k] == length + large);
    1042                     next[kRow] = k;
    1043                     previous[k] = kRow;
    1044                }
    1045           }
    1046      }
    1047      // do rest
    1048      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1049           if (permuteInverse_[iRow] < 0)
    1050                permuteInverse_[iRow] = done++;
    1051      }
    1052      COIN_DETAIL_PRINT(printf("%d compressions, %d expansions\n",
    1053                               numberCompressions, numberExpansions));
    1054      assert (done == numberRows_);
    1055      delete [] sort;
    1056      delete [] order;
    1057      delete [] which;
    1058      delete [] mark;
    1059      delete [] first;
    1060      delete [] next;
    1061      delete [] previous;
    1062      delete [] choleskyRow_;
    1063      choleskyRow_ = NULL;
    1064      delete [] choleskyStart_;
    1065      choleskyStart_ = NULL;
     787  permuteInverse_ = new int[numberRows_];
     788  permute_ = new int[numberRows_];
     789  // Do ordering
     790  int returnCode = 0;
     791  // get more space and full matrix
     792  int space = 6 * sizeFactor_ + 100000;
     793  int *temp = new int[space];
     794  int *which = new int[2 * numberRows_];
     795  CoinBigIndex *tempStart = new CoinBigIndex[numberRows_ + 1];
     796  memset(which, 0, numberRows_ * sizeof(int));
     797  for (int iRow = 0; iRow < numberRows_; iRow++) {
     798    which[iRow] += choleskyStart_[iRow + 1] - choleskyStart_[iRow] - 1;
     799    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
     800      int jRow = choleskyRow_[j];
     801      which[jRow]++;
     802    }
     803  }
     804  CoinBigIndex sizeFactor = 0;
     805  for (int iRow = 0; iRow < numberRows_; iRow++) {
     806    int length = which[iRow];
     807    permute_[iRow] = length;
     808    tempStart[iRow] = sizeFactor;
     809    which[iRow] = sizeFactor;
     810    sizeFactor += length;
     811  }
     812  for (int iRow = 0; iRow < numberRows_; iRow++) {
     813    assert(choleskyRow_[choleskyStart_[iRow]] == iRow);
     814    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
     815      int jRow = choleskyRow_[j];
     816      int put = which[iRow];
     817      temp[put] = jRow;
     818      which[iRow]++;
     819      put = which[jRow];
     820      temp[put] = iRow;
     821      which[jRow]++;
     822    }
     823  }
     824  for (int iRow = 1; iRow < numberRows_; iRow++)
     825    assert(which[iRow - 1] == tempStart[iRow]);
     826  CoinBigIndex lastSpace = sizeFactor;
     827  delete[] choleskyRow_;
     828  choleskyRow_ = temp;
     829  delete[] choleskyStart_;
     830  choleskyStart_ = tempStart;
     831  // linked lists of sizes and lengths
     832  int *first = new int[numberRows_];
     833  int *next = new int[numberRows_];
     834  int *previous = new int[numberRows_];
     835  char *mark = new char[numberRows_];
     836  memset(mark, 0, numberRows_);
     837  CoinBigIndex *sort = new CoinBigIndex[numberRows_];
     838  int *order = new int[numberRows_];
     839  for (int iRow = 0; iRow < numberRows_; iRow++) {
     840    first[iRow] = -1;
     841    next[iRow] = -1;
     842    previous[iRow] = -1;
     843    permuteInverse_[iRow] = -1;
     844  }
     845  int large = 1000 + 2 * numberRows_;
     846  for (int iRow = 0; iRow < numberRows_; iRow++) {
     847    int n = permute_[iRow];
     848    if (first[n] < 0) {
     849      first[n] = iRow;
     850      previous[iRow] = n + large;
     851      next[iRow] = n + 2 * large;
     852    } else {
     853      int k = first[n];
     854      assert(k < numberRows_);
     855      first[n] = iRow;
     856      previous[iRow] = n + large;
     857      assert(previous[k] == n + large);
     858      next[iRow] = k;
     859      previous[k] = iRow;
     860    }
     861  }
     862  int smallest = 0;
     863  int done = 0;
     864  int numberCompressions = 0;
     865  int numberExpansions = 0;
     866  while (smallest < numberRows_) {
     867    if (first[smallest] < 0 || first[smallest] > numberRows_) {
     868      smallest++;
     869      continue;
     870    }
     871    int iRow = first[smallest];
     872    if (false) {
     873      int kRow = -1;
     874      int ss = 999999;
     875      for (int jRow = numberRows_ - 1; jRow >= 0; jRow--) {
     876        if (permuteInverse_[jRow] < 0) {
     877          int length = permute_[jRow];
     878          assert(length > 0);
     879          for (CoinBigIndex j = choleskyStart_[jRow];
     880               j < choleskyStart_[jRow] + length; j++) {
     881            int jjRow = choleskyRow_[j];
     882            assert(permuteInverse_[jjRow] < 0);
     883          }
     884          if (length < ss) {
     885            ss = length;
     886            kRow = jRow;
     887          }
     888        }
     889      }
     890      assert(smallest == ss);
     891      printf("list chose %d - full chose %d - length %d\n",
     892        iRow, kRow, ss);
     893    }
     894    int kNext = next[iRow];
     895    first[smallest] = kNext;
     896    if (kNext < numberRows_)
     897      previous[kNext] = previous[iRow];
     898    previous[iRow] = -1;
     899    next[iRow] = -1;
     900    permuteInverse_[iRow] = done;
     901    done++;
     902    // Now add edges
     903    CoinBigIndex start = choleskyStart_[iRow];
     904    CoinBigIndex end = choleskyStart_[iRow] + permute_[iRow];
     905    int nSave = 0;
     906    for (CoinBigIndex k = start; k < end; k++) {
     907      int kRow = choleskyRow_[k];
     908      assert(permuteInverse_[kRow] < 0);
     909      //if (permuteInverse_[kRow]<0)
     910      which[nSave++] = kRow;
     911    }
     912    for (int i = 0; i < nSave; i++) {
     913      int kRow = which[i];
     914      int length = permute_[kRow];
     915      CoinBigIndex start = choleskyStart_[kRow];
     916      CoinBigIndex end = choleskyStart_[kRow] + length;
     917      for (CoinBigIndex j = start; j < end; j++) {
     918        int jRow = choleskyRow_[j];
     919        mark[jRow] = 1;
     920      }
     921      mark[kRow] = 1;
     922      int n = nSave;
     923      for (int j = 0; j < nSave; j++) {
     924        int jRow = which[j];
     925        if (!mark[jRow]) {
     926          which[n++] = jRow;
     927        }
     928      }
     929      for (CoinBigIndex j = start; j < end; j++) {
     930        int jRow = choleskyRow_[j];
     931        mark[jRow] = 0;
     932      }
     933      mark[kRow] = 0;
     934      if (n > nSave) {
     935        bool inPlace = (n - nSave == 1);
     936        if (!inPlace) {
     937          // extend
     938          int length = n - nSave + end - start;
     939          if (length + lastSpace > space) {
     940            // need to compress
     941            numberCompressions++;
     942            int newN = 0;
     943            for (int iRow = 0; iRow < numberRows_; iRow++) {
     944              CoinBigIndex start = choleskyStart_[iRow];
     945              if (permuteInverse_[iRow] < 0) {
     946                sort[newN] = start;
     947                order[newN++] = iRow;
     948              } else {
     949                choleskyStart_[iRow] = -1;
     950                permute_[iRow] = 0;
     951              }
     952            }
     953            CoinSort_2(sort, sort + newN, order);
     954            sizeFactor = 0;
     955            for (int k = 0; k < newN; k++) {
     956              int iRow = order[k];
     957              int length = permute_[iRow];
     958              CoinBigIndex start = choleskyStart_[iRow];
     959              choleskyStart_[iRow] = sizeFactor;
     960              for (int j = 0; j < length; j++)
     961                choleskyRow_[sizeFactor + j] = choleskyRow_[start + j];
     962              sizeFactor += length;
     963            }
     964            lastSpace = sizeFactor;
     965            if ((sizeFactor + numberRows_ + 1000) * 4 > 3 * space) {
     966              // need to expand
     967              numberExpansions++;
     968              space = (3 * space) / 2;
     969              int *temp = new int[space];
     970              memcpy(temp, choleskyRow_, sizeFactor * sizeof(int));
     971              delete[] choleskyRow_;
     972              choleskyRow_ = temp;
     973            }
     974          }
     975        }
     976        // Now add
     977        start = choleskyStart_[kRow];
     978        end = choleskyStart_[kRow] + permute_[kRow];
     979        if (!inPlace)
     980          choleskyStart_[kRow] = lastSpace;
     981        CoinBigIndex put = choleskyStart_[kRow];
     982        for (CoinBigIndex j = start; j < end; j++) {
     983          int jRow = choleskyRow_[j];
     984          if (permuteInverse_[jRow] < 0)
     985            choleskyRow_[put++] = jRow;
     986        }
     987        for (int j = nSave; j < n; j++) {
     988          int jRow = which[j];
     989          choleskyRow_[put++] = jRow;
     990        }
     991        if (!inPlace) {
     992          permute_[kRow] = put - lastSpace;
     993          lastSpace = put;
     994        } else {
     995          permute_[kRow] = put - choleskyStart_[kRow];
     996        }
     997      } else {
     998        // pack down for new counts
     999        CoinBigIndex put = start;
     1000        for (CoinBigIndex j = start; j < end; j++) {
     1001          int jRow = choleskyRow_[j];
     1002          if (permuteInverse_[jRow] < 0)
     1003            choleskyRow_[put++] = jRow;
     1004        }
     1005        permute_[kRow] = put - start;
     1006      }
     1007      // take out
     1008      int iNext = next[kRow];
     1009      int iPrevious = previous[kRow];
     1010      if (iPrevious < numberRows_) {
     1011        next[iPrevious] = iNext;
     1012      } else {
     1013        assert(iPrevious == length + large);
     1014        first[length] = iNext;
     1015      }
     1016      if (iNext < numberRows_) {
     1017        previous[iNext] = iPrevious;
     1018      } else {
     1019        assert(iNext == length + 2 * large);
     1020      }
     1021      // put in
     1022      length = permute_[kRow];
     1023      smallest = CoinMin(smallest, length);
     1024      if (first[length] < 0 || first[length] > numberRows_) {
     1025        first[length] = kRow;
     1026        previous[kRow] = length + large;
     1027        next[kRow] = length + 2 * large;
     1028      } else {
     1029        int k = first[length];
     1030        assert(k < numberRows_);
     1031        first[length] = kRow;
     1032        previous[kRow] = length + large;
     1033        assert(previous[k] == length + large);
     1034        next[kRow] = k;
     1035        previous[k] = kRow;
     1036      }
     1037    }
     1038  }
     1039  // do rest
     1040  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1041    if (permuteInverse_[iRow] < 0)
     1042      permuteInverse_[iRow] = done++;
     1043  }
     1044  COIN_DETAIL_PRINT(printf("%d compressions, %d expansions\n",
     1045    numberCompressions, numberExpansions));
     1046  assert(done == numberRows_);
     1047  delete[] sort;
     1048  delete[] order;
     1049  delete[] which;
     1050  delete[] mark;
     1051  delete[] first;
     1052  delete[] next;
     1053  delete[] previous;
     1054  delete[] choleskyRow_;
     1055  choleskyRow_ = NULL;
     1056  delete[] choleskyStart_;
     1057  choleskyStart_ = NULL;
    10661058#ifndef NDEBUG
    1067      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1068           permute_[iRow] = -1;
    1069      }
    1070 #endif
    1071      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1072           permute_[permuteInverse_[iRow]] = iRow;
    1073      }
     1059  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1060    permute_[iRow] = -1;
     1061  }
     1062#endif
     1063  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1064    permute_[permuteInverse_[iRow]] = iRow;
     1065  }
    10741066#ifndef NDEBUG
    1075      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1076           assert(permute_[iRow] >= 0 && permute_[iRow] < numberRows_);
    1077      }
    1078 #endif
    1079      return returnCode;
     1067  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1068    assert(permute_[iRow] >= 0 && permute_[iRow] < numberRows_);
     1069  }
     1070#endif
     1071  return returnCode;
    10801072}
    1081 #elif BASE_ORDER==2
     1073#elif BASE_ORDER == 2
    10821074/*----------------------------------------------------------------------------*/
    10831075/*      (C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
     
    10981090#include <stdlib.h>
    10991091
    1100 typedef int             WSI;
     1092typedef int WSI;
    11011093
    1102 #define NORDTHRESH      7
    1103 #define MAXIW           2147000000
     1094#define NORDTHRESH 7
     1095#define MAXIW 2147000000
    11041096
    11051097//#define WSSMP_DBG
    11061098#ifdef WSSMP_DBG
    1107 static void chk (WSI ind, WSI i, WSI lim)
     1099static void chk(WSI ind, WSI i, WSI lim)
    11081100{
    11091101
    1110      if (i <= 0 || i > lim) {
    1111           printf("%d: chk: myamlf: WAR**: i, lim = %d, %d\n", ind, i, lim);
    1112      }
     1102  if (i <= 0 || i > lim) {
     1103    printf("%d: chk: myamlf: WAR**: i, lim = %d, %d\n", ind, i, lim);
     1104  }
    11131105}
    11141106#endif
     
    11161108static void
    11171109myamlf(WSI n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI varbl[],
    1118        WSI snxt[], WSI perm[], WSI invp[], WSI head[], WSI lsize[],
    1119        WSI flag[], WSI erscore[], WSI locaux, WSI adjln, WSI speed)
     1110  WSI snxt[], WSI perm[], WSI invp[], WSI head[], WSI lsize[],
     1111  WSI flag[], WSI erscore[], WSI locaux, WSI adjln, WSI speed)
    11201112{
    1121      WSI l, i, j, k;
    1122      double x, y;
    1123      WSI maxmum, fltag, nodeg, scln, nm1, deg, tn,
    1124          locatns, ipp, jpp, nnode, numpiv, node,
    1125          nodeln, currloc, counter, numii, mindeg,
    1126          i0, i1, i2, i4, i5, i6, i7, i9,
    1127          j0, j1, j2, j3, j4, j5, j6, j7, j8, j9;
    1128      /*                                             n cannot be less than NORDTHRESH
     1113  WSI l, i, j, k;
     1114  double x, y;
     1115  WSI maxmum, fltag, nodeg, scln, nm1, deg, tn,
     1116    locatns, ipp, jpp, nnode, numpiv, node,
     1117    nodeln, currloc, counter, numii, mindeg,
     1118    i0, i1, i2, i4, i5, i6, i7, i9,
     1119    j0, j1, j2, j3, j4, j5, j6, j7, j8, j9;
     1120  /*                                             n cannot be less than NORDTHRESH
    11291121      if (n < 3) {
    11301122         if (n > 0) perm[0] = invp[0] = 1;
     
    11341126     */
    11351127#ifdef WSSMP_DBG
    1136      printf("Myamlf: n, locaux, adjln, speed = %d,%d,%d,%d\n", n, locaux, adjln, speed);
    1137      for (i = 0; i < n; i++) flag[i] = 0;
    1138      k = xadj[0];
    1139      for (i = 1; i <= n; i++) {
    1140           j = k + dgree[i-1];
    1141           if (j < k || k < 1) printf("ERR**: myamlf: i, j, k = %d,%d,%d\n", i, j, k);
    1142           for (l = k; l < j; l++) {
    1143                if (adjncy[l - 1] < 1 || adjncy[l - 1] > n || adjncy[l - 1] == i)
    1144                     printf("ERR**: myamlf: i, l, adjj[l] = %d,%d,%d\n", i, l, adjncy[l - 1]);
    1145                if (flag[adjncy[l - 1] - 1] == i)
    1146                     printf("ERR**: myamlf: duplicate entry: %d - %d\n", i, adjncy[l - 1]);
    1147                flag[adjncy[l - 1] - 1] = i;
    1148                nm1 = adjncy[l - 1] - 1;
    1149                for (fltag = xadj[nm1]; fltag < xadj[nm1] + dgree[nm1]; fltag++) {
    1150                     if (adjncy[fltag - 1] == i) goto L100;
    1151                }
    1152                printf("ERR**: Unsymmetric %d %d\n", i, fltag);
    1153 L100:
    1154                ;
    1155           }
    1156           k = j;
    1157           if (k > locaux) printf("ERR**: i, j, k, locaux = %d, %d, %d, %d\n",
    1158                                       i, j, k, locaux);
    1159      }
    1160      if (n*(n - 1) < locaux - 1) printf("ERR**: myamlf: too many nnz, n, ne = %d, %d\n",
    1161                                              n, adjln);
    1162      for (i = 0; i < n; i++) flag[i] = 1;
    1163      if (n > 10000) printf ("Finished error checking in AMF\n");
    1164      for (i = locaux; i <= adjln; i++) adjncy[i - 1] = -22;
     1128  printf("Myamlf: n, locaux, adjln, speed = %d,%d,%d,%d\n", n, locaux, adjln, speed);
     1129  for (i = 0; i < n; i++)
     1130    flag[i] = 0;
     1131  k = xadj[0];
     1132  for (i = 1; i <= n; i++) {
     1133    j = k + dgree[i - 1];
     1134    if (j < k || k < 1)
     1135      printf("ERR**: myamlf: i, j, k = %d,%d,%d\n", i, j, k);
     1136    for (l = k; l < j; l++) {
     1137      if (adjncy[l - 1] < 1 || adjncy[l - 1] > n || adjncy[l - 1] == i)
     1138        printf("ERR**: myamlf: i, l, adjj[l] = %d,%d,%d\n", i, l, adjncy[l - 1]);
     1139      if (flag[adjncy[l - 1] - 1] == i)
     1140        printf("ERR**: myamlf: duplicate entry: %d - %d\n", i, adjncy[l - 1]);
     1141      flag[adjncy[l - 1] - 1] = i;
     1142      nm1 = adjncy[l - 1] - 1;
     1143      for (fltag = xadj[nm1]; fltag < xadj[nm1] + dgree[nm1]; fltag++) {
     1144        if (adjncy[fltag - 1] == i)
     1145          goto L100;
     1146      }
     1147      printf("ERR**: Unsymmetric %d %d\n", i, fltag);
     1148    L100:;
     1149    }
     1150    k = j;
     1151    if (k > locaux)
     1152      printf("ERR**: i, j, k, locaux = %d, %d, %d, %d\n",
     1153        i, j, k, locaux);
     1154  }
     1155  if (n * (n - 1) < locaux - 1)
     1156    printf("ERR**: myamlf: too many nnz, n, ne = %d, %d\n",
     1157      n, adjln);
     1158  for (i = 0; i < n; i++)
     1159    flag[i] = 1;
     1160  if (n > 10000)
     1161    printf("Finished error checking in AMF\n");
     1162  for (i = locaux; i <= adjln; i++)
     1163    adjncy[i - 1] = -22;
    11651164#endif
    11661165
    1167      maxmum = 0;
    1168      mindeg = 1;
    1169      fltag = 2;
    1170      locatns = locaux - 1;
    1171      nm1 = n - 1;
    1172      counter = 1;
    1173      for (l = 0; l < n; l++) {
    1174           j = erscore[l];
    1175 #ifdef WSSMP_DBG
    1176           chk(1, j, n);
    1177 #endif
    1178           if (j > 0) {
    1179                nnode = head[j - 1];
    1180                if (nnode) perm[nnode - 1] = l + 1;
    1181                snxt[l] = nnode;
    1182                head[j - 1] = l + 1;
     1166  maxmum = 0;
     1167  mindeg = 1;
     1168  fltag = 2;
     1169  locatns = locaux - 1;
     1170  nm1 = n - 1;
     1171  counter = 1;
     1172  for (l = 0; l < n; l++) {
     1173    j = erscore[l];
     1174#ifdef WSSMP_DBG
     1175    chk(1, j, n);
     1176#endif
     1177    if (j > 0) {
     1178      nnode = head[j - 1];
     1179      if (nnode)
     1180        perm[nnode - 1] = l + 1;
     1181      snxt[l] = nnode;
     1182      head[j - 1] = l + 1;
     1183    } else {
     1184      invp[l] = -(counter++);
     1185      flag[l] = xadj[l] = 0;
     1186    }
     1187  }
     1188  while (counter <= n) {
     1189    for (deg = mindeg; head[deg - 1] < 1; deg++) {
     1190    };
     1191    nodeg = 0;
     1192#ifdef WSSMP_DBG
     1193    chk(2, deg, n - 1);
     1194#endif
     1195    node = head[deg - 1];
     1196#ifdef WSSMP_DBG
     1197    chk(3, node, n);
     1198#endif
     1199    mindeg = deg;
     1200    nnode = snxt[node - 1];
     1201    if (nnode)
     1202      perm[nnode - 1] = 0;
     1203    head[deg - 1] = nnode;
     1204    nodeln = invp[node - 1];
     1205    numpiv = varbl[node - 1];
     1206    invp[node - 1] = -counter;
     1207    counter += numpiv;
     1208    varbl[node - 1] = -numpiv;
     1209    if (nodeln) {
     1210      j4 = locaux;
     1211      i5 = lsize[node - 1] - nodeln;
     1212      i2 = nodeln + 1;
     1213      l = xadj[node - 1];
     1214      for (i6 = 1; i6 <= i2; ++i6) {
     1215        if (i6 == i2) {
     1216          tn = node, i0 = l, scln = i5;
     1217        } else {
     1218#ifdef WSSMP_DBG
     1219          chk(4, l, adjln);
     1220#endif
     1221          tn = adjncy[l - 1];
     1222          l++;
     1223#ifdef WSSMP_DBG
     1224          chk(5, tn, n);
     1225#endif
     1226          i0 = xadj[tn - 1];
     1227          scln = lsize[tn - 1];
     1228        }
     1229        for (i7 = 1; i7 <= scln; ++i7) {
     1230#ifdef WSSMP_DBG
     1231          chk(6, i0, adjln);
     1232#endif
     1233          i = adjncy[i0 - 1];
     1234          i0++;
     1235#ifdef WSSMP_DBG
     1236          chk(7, i, n);
     1237#endif
     1238          numii = varbl[i - 1];
     1239          if (numii > 0) {
     1240            if (locaux > adjln) {
     1241              lsize[node - 1] -= i6;
     1242              xadj[node - 1] = l;
     1243              if (!lsize[node - 1])
     1244                xadj[node - 1] = 0;
     1245              xadj[tn - 1] = i0;
     1246              lsize[tn - 1] = scln - i7;
     1247              if (!lsize[tn - 1])
     1248                xadj[tn - 1] = 0;
     1249              for (j = 1; j <= n; j++) {
     1250                i4 = xadj[j - 1];
     1251                if (i4 > 0) {
     1252                  xadj[j - 1] = adjncy[i4 - 1];
     1253#ifdef WSSMP_DBG
     1254                  chk(8, i4, adjln);
     1255#endif
     1256                  adjncy[i4 - 1] = -j;
     1257                }
     1258              }
     1259              i9 = j4 - 1;
     1260              j6 = j7 = 1;
     1261              while (j6 <= i9) {
     1262#ifdef WSSMP_DBG
     1263                chk(9, j6, adjln);
     1264#endif
     1265                j = -adjncy[j6 - 1];
     1266                j6++;
     1267                if (j > 0) {
     1268#ifdef WSSMP_DBG
     1269                  chk(10, j7, adjln);
     1270                  chk(11, j, n);
     1271#endif
     1272                  adjncy[j7 - 1] = xadj[j - 1];
     1273                  xadj[j - 1] = j7++;
     1274                  j8 = lsize[j - 1] - 1 + j7;
     1275                  for (; j7 < j8; j7++, j6++)
     1276                    adjncy[j7 - 1] = adjncy[j6 - 1];
     1277                }
     1278              }
     1279              for (j0 = j7; j4 < locaux; j4++, j7++) {
     1280#ifdef WSSMP_DBG
     1281                chk(12, j4, adjln);
     1282#endif
     1283                adjncy[j7 - 1] = adjncy[j4 - 1];
     1284              }
     1285              j4 = j0;
     1286              locaux = j7;
     1287              i0 = xadj[tn - 1];
     1288              l = xadj[node - 1];
     1289            }
     1290            adjncy[locaux - 1] = i;
     1291            locaux++;
     1292            varbl[i - 1] = -numii;
     1293            nodeg += numii;
     1294            ipp = perm[i - 1];
     1295            nnode = snxt[i - 1];
     1296#ifdef WSSMP_DBG
     1297            if (ipp)
     1298              chk(13, ipp, n);
     1299            if (nnode)
     1300              chk(14, nnode, n);
     1301#endif
     1302            if (ipp)
     1303              snxt[ipp - 1] = nnode;
     1304            else
     1305              head[erscore[i - 1] - 1] = nnode;
     1306            if (nnode)
     1307              perm[nnode - 1] = ipp;
     1308          }
     1309        }
     1310        if (tn != node) {
     1311          flag[tn - 1] = 0, xadj[tn - 1] = -node;
     1312        }
     1313      }
     1314      currloc = (j5 = locaux) - j4;
     1315      locatns += currloc;
     1316    } else {
     1317      i1 = (j4 = xadj[node - 1]) + lsize[node - 1];
     1318      for (j = j5 = j4; j < i1; j++) {
     1319#ifdef WSSMP_DBG
     1320        chk(15, j, adjln);
     1321#endif
     1322        i = adjncy[j - 1];
     1323#ifdef WSSMP_DBG
     1324        chk(16, i, n);
     1325#endif
     1326        numii = varbl[i - 1];
     1327        if (numii > 0) {
     1328          nodeg += numii;
     1329          varbl[i - 1] = -numii;
     1330#ifdef WSSMP_DBG
     1331          chk(17, j5, adjln);
     1332#endif
     1333          adjncy[j5 - 1] = i;
     1334          ipp = perm[i - 1];
     1335          nnode = snxt[i - 1];
     1336          j5++;
     1337#ifdef WSSMP_DBG
     1338          if (ipp)
     1339            chk(18, ipp, n);
     1340          if (nnode)
     1341            chk(19, nnode, n);
     1342#endif
     1343          if (ipp)
     1344            snxt[ipp - 1] = nnode;
     1345          else
     1346            head[erscore[i - 1] - 1] = nnode;
     1347          if (nnode)
     1348            perm[nnode - 1] = ipp;
     1349        }
     1350      }
     1351      currloc = 0;
     1352    }
     1353#ifdef WSSMP_DBG
     1354    chk(20, node, n);
     1355#endif
     1356    lsize[node - 1] = j5 - (xadj[node - 1] = j4);
     1357    dgree[node - 1] = nodeg;
     1358    if (fltag + n < 0 || fltag + n > MAXIW) {
     1359      for (i = 1; i <= n; i++)
     1360        if (flag[i - 1])
     1361          flag[i - 1] = 1;
     1362      fltag = 2;
     1363    }
     1364    for (j3 = j4; j3 < j5; j3++) {
     1365#ifdef WSSMP_DBG
     1366      chk(21, j3, adjln);
     1367#endif
     1368      i = adjncy[j3 - 1];
     1369#ifdef WSSMP_DBG
     1370      chk(22, i, n);
     1371#endif
     1372      j = invp[i - 1];
     1373      if (j > 0) {
     1374        i4 = fltag - (numii = -varbl[i - 1]);
     1375        i2 = xadj[i - 1] + j;
     1376        for (l = xadj[i - 1]; l < i2; l++) {
     1377#ifdef WSSMP_DBG
     1378          chk(23, l, adjln);
     1379#endif
     1380          tn = adjncy[l - 1];
     1381#ifdef WSSMP_DBG
     1382          chk(24, tn, n);
     1383#endif
     1384          j9 = flag[tn - 1];
     1385          if (j9 >= fltag)
     1386            j9 -= numii;
     1387          else if (j9)
     1388            j9 = dgree[tn - 1] + i4;
     1389          flag[tn - 1] = j9;
     1390        }
     1391      }
     1392    }
     1393    for (j3 = j4; j3 < j5; j3++) {
     1394#ifdef WSSMP_DBG
     1395      chk(25, j3, adjln);
     1396#endif
     1397      i = adjncy[j3 - 1];
     1398      i5 = deg = 0;
     1399#ifdef WSSMP_DBG
     1400      chk(26, i, n);
     1401#endif
     1402      j1 = (i4 = xadj[i - 1]) + invp[i - 1];
     1403      for (l = j0 = i4; l < j1; l++) {
     1404#ifdef WSSMP_DBG
     1405        chk(27, l, adjln);
     1406#endif
     1407        tn = adjncy[l - 1];
     1408#ifdef WSSMP_DBG
     1409        chk(70, tn, n);
     1410#endif
     1411        j8 = flag[tn - 1];
     1412        if (j8) {
     1413          deg += (j8 - fltag);
     1414#ifdef WSSMP_DBG
     1415          chk(28, i4, adjln);
     1416#endif
     1417          adjncy[i4 - 1] = tn;
     1418          i5 += tn;
     1419          i4++;
     1420          while (i5 >= nm1)
     1421            i5 -= nm1;
     1422        }
     1423      }
     1424      invp[i - 1] = (j2 = i4) - j0 + 1;
     1425      i2 = j0 + lsize[i - 1];
     1426      for (l = j1; l < i2; l++) {
     1427#ifdef WSSMP_DBG
     1428        chk(29, l, adjln);
     1429#endif
     1430        j = adjncy[l - 1];
     1431#ifdef WSSMP_DBG
     1432        chk(30, j, n);
     1433#endif
     1434        numii = varbl[j - 1];
     1435        if (numii > 0) {
     1436          deg += numii;
     1437#ifdef WSSMP_DBG
     1438          chk(31, i4, adjln);
     1439#endif
     1440          adjncy[i4 - 1] = j;
     1441          i5 += j;
     1442          i4++;
     1443          while (i5 >= nm1)
     1444            i5 -= nm1;
     1445        }
     1446      }
     1447      if (invp[i - 1] == 1 && j2 == i4) {
     1448        numii = -varbl[i - 1];
     1449        xadj[i - 1] = -node;
     1450        nodeg -= numii;
     1451        counter += numii;
     1452        numpiv += numii;
     1453        invp[i - 1] = varbl[i - 1] = 0;
     1454      } else {
     1455        if (dgree[i - 1] > deg)
     1456          dgree[i - 1] = deg;
     1457#ifdef WSSMP_DBG
     1458        chk(32, j2, adjln);
     1459        chk(33, j0, adjln);
     1460#endif
     1461        adjncy[i4 - 1] = adjncy[j2 - 1];
     1462        adjncy[j2 - 1] = adjncy[j0 - 1];
     1463        adjncy[j0 - 1] = node;
     1464        lsize[i - 1] = i4 - j0 + 1;
     1465        i5++;
     1466#ifdef WSSMP_DBG
     1467        chk(35, i5, n);
     1468#endif
     1469        j = head[i5 - 1];
     1470        if (j > 0) {
     1471          snxt[i - 1] = perm[j - 1];
     1472          perm[j - 1] = i;
     1473        } else {
     1474          snxt[i - 1] = -j;
     1475          head[i5 - 1] = -i;
     1476        }
     1477        perm[i - 1] = i5;
     1478      }
     1479    }
     1480#ifdef WSSMP_DBG
     1481    chk(36, node, n);
     1482#endif
     1483    dgree[node - 1] = nodeg;
     1484    if (maxmum < nodeg)
     1485      maxmum = nodeg;
     1486    fltag += maxmum;
     1487#ifdef WSSMP_DBG
     1488    if (fltag + n + 1 < 0)
     1489      printf("ERR**: fltag = %d (A)\n", fltag);
     1490#endif
     1491    for (j3 = j4; j3 < j5; ++j3) {
     1492#ifdef WSSMP_DBG
     1493      chk(37, j3, adjln);
     1494#endif
     1495      i = adjncy[j3 - 1];
     1496#ifdef WSSMP_DBG
     1497      chk(38, i, n);
     1498#endif
     1499      if (varbl[i - 1] < 0) {
     1500        i5 = perm[i - 1];
     1501#ifdef WSSMP_DBG
     1502        chk(39, i5, n);
     1503#endif
     1504        j = head[i5 - 1];
     1505        if (j) {
     1506          if (j < 0) {
     1507            head[i5 - 1] = 0, i = -j;
    11831508          } else {
    1184                invp[l] = -(counter++);
    1185                flag[l] = xadj[l] = 0;
    1186           }
    1187      }
    1188      while (counter <= n) {
    1189           for (deg = mindeg; head[deg - 1] < 1; deg++) {};
    1190           nodeg = 0;
    1191 #ifdef WSSMP_DBG
    1192           chk(2, deg, n - 1);
    1193 #endif
    1194           node = head[deg - 1];
    1195 #ifdef WSSMP_DBG
    1196           chk(3, node, n);
    1197 #endif
    1198           mindeg = deg;
    1199           nnode = snxt[node - 1];
    1200           if (nnode) perm[nnode - 1] = 0;
    1201           head[deg - 1] = nnode;
    1202           nodeln = invp[node - 1];
    1203           numpiv = varbl[node - 1];
    1204           invp[node - 1] = -counter;
    1205           counter += numpiv;
    1206           varbl[node - 1] = -numpiv;
    1207           if (nodeln) {
    1208                j4 = locaux;
    1209                i5 = lsize[node - 1] - nodeln;
    1210                i2 = nodeln + 1;
    1211                l = xadj[node - 1];
    1212                for (i6 = 1; i6 <= i2; ++i6) {
    1213                     if (i6 == i2) {
    1214                          tn = node, i0 = l, scln = i5;
    1215                     } else {
    1216 #ifdef WSSMP_DBG
    1217                          chk(4, l, adjln);
    1218 #endif
    1219                          tn = adjncy[l-1];
    1220                          l++;
    1221 #ifdef WSSMP_DBG
    1222                          chk(5, tn, n);
    1223 #endif
    1224                          i0 = xadj[tn - 1];
    1225                          scln = lsize[tn - 1];
     1509#ifdef WSSMP_DBG
     1510            chk(40, j, n);
     1511#endif
     1512            i = perm[j - 1];
     1513            perm[j - 1] = 0;
     1514          }
     1515          while (i) {
     1516#ifdef WSSMP_DBG
     1517            chk(41, i, n);
     1518#endif
     1519            if (!snxt[i - 1]) {
     1520              i = 0;
     1521            } else {
     1522              k = invp[i - 1];
     1523              i2 = xadj[i - 1] + (scln = lsize[i - 1]);
     1524              for (l = xadj[i - 1] + 1; l < i2; l++) {
     1525#ifdef WSSMP_DBG
     1526                chk(42, l, adjln);
     1527                chk(43, adjncy[l - 1], n);
     1528#endif
     1529                flag[adjncy[l - 1] - 1] = fltag;
     1530              }
     1531              jpp = i;
     1532              j = snxt[i - 1];
     1533              while (j) {
     1534#ifdef WSSMP_DBG
     1535                chk(44, j, n);
     1536#endif
     1537                if (lsize[j - 1] == scln && invp[j - 1] == k) {
     1538                  i2 = xadj[j - 1] + scln;
     1539                  j8 = 1;
     1540                  for (l = xadj[j - 1] + 1; l < i2; l++) {
     1541#ifdef WSSMP_DBG
     1542                    chk(45, l, adjln);
     1543                    chk(46, adjncy[l - 1], n);
     1544#endif
     1545                    if (flag[adjncy[l - 1] - 1] != fltag) {
     1546                      j8 = 0;
     1547                      break;
    12261548                    }
    1227                     for (i7 = 1; i7 <= scln; ++i7) {
    1228 #ifdef WSSMP_DBG
    1229                          chk(6, i0, adjln);
    1230 #endif
    1231                          i = adjncy[i0 - 1];
    1232                          i0++;
    1233 #ifdef WSSMP_DBG
    1234                          chk(7, i, n);
    1235 #endif
    1236                          numii = varbl[i - 1];
    1237                          if (numii > 0) {
    1238                               if (locaux > adjln) {
    1239                                    lsize[node - 1] -= i6;
    1240                                    xadj[node - 1] = l;
    1241                                    if (!lsize[node - 1]) xadj[node - 1] = 0;
    1242                                    xadj[tn - 1] = i0;
    1243                                    lsize[tn - 1] = scln - i7;
    1244                                    if (!lsize[tn - 1]) xadj[tn - 1] = 0;
    1245                                    for (j = 1; j <= n; j++) {
    1246                                         i4 = xadj[j - 1];
    1247                                         if (i4 > 0) {
    1248                                              xadj[j - 1] = adjncy[i4 - 1];
    1249 #ifdef WSSMP_DBG
    1250                                              chk(8, i4, adjln);
    1251 #endif
    1252                                              adjncy[i4 - 1] = -j;
    1253                                         }
    1254                                    }
    1255                                    i9 = j4 - 1;
    1256                                    j6 = j7 = 1;
    1257                                    while (j6 <= i9) {
    1258 #ifdef WSSMP_DBG
    1259                                         chk(9, j6, adjln);
    1260 #endif
    1261                                         j = -adjncy[j6 - 1];
    1262                                         j6++;
    1263                                         if (j > 0) {
    1264 #ifdef WSSMP_DBG
    1265                                              chk(10, j7, adjln);
    1266                                              chk(11, j, n);
    1267 #endif
    1268                                              adjncy[j7 - 1] = xadj[j - 1];
    1269                                              xadj[j - 1] = j7++;
    1270                                              j8 = lsize[j - 1] - 1 + j7;
    1271                                              for (; j7 < j8; j7++, j6++) adjncy[j7-1] = adjncy[j6-1];
    1272                                         }
    1273                                    }
    1274                                    for (j0 = j7; j4 < locaux; j4++, j7++) {
    1275 #ifdef WSSMP_DBG
    1276                                         chk(12, j4, adjln);
    1277 #endif
    1278                                         adjncy[j7 - 1] = adjncy[j4 - 1];
    1279                                    }
    1280                                    j4 = j0;
    1281                                    locaux = j7;
    1282                                    i0 = xadj[tn - 1];
    1283                                    l = xadj[node - 1];
    1284                               }
    1285                               adjncy[locaux-1] = i;
    1286                               locaux++;
    1287                               varbl[i - 1] = -numii;
    1288                               nodeg += numii;
    1289                               ipp = perm[i - 1];
    1290                               nnode = snxt[i - 1];
    1291 #ifdef WSSMP_DBG
    1292                               if (ipp) chk(13, ipp, n);
    1293                               if (nnode) chk(14, nnode, n);
    1294 #endif
    1295                               if (ipp) snxt[ipp - 1] = nnode;
    1296                               else head[erscore[i - 1] - 1] = nnode;
    1297                               if (nnode) perm[nnode - 1] = ipp;
    1298                          }
    1299                     }
    1300                     if (tn != node) {
    1301                          flag[tn - 1] = 0, xadj[tn - 1] = -node;
    1302                     }
    1303                }
    1304                currloc = (j5 = locaux) - j4;
    1305                locatns += currloc;
    1306           } else {
    1307                i1 = (j4 = xadj[node - 1]) + lsize[node - 1];
    1308                for (j = j5 = j4; j < i1; j++) {
    1309 #ifdef WSSMP_DBG
    1310                     chk(15, j, adjln);
    1311 #endif
    1312                     i = adjncy[j - 1];
    1313 #ifdef WSSMP_DBG
    1314                     chk(16, i, n);
    1315 #endif
    1316                     numii = varbl[i - 1];
    1317                     if (numii > 0) {
    1318                          nodeg += numii;
    1319                          varbl[i - 1] = -numii;
    1320 #ifdef WSSMP_DBG
    1321                          chk(17, j5, adjln);
    1322 #endif
    1323                          adjncy[j5-1] = i;
    1324                          ipp = perm[i - 1];
    1325                          nnode = snxt[i - 1];
    1326                          j5++;
    1327 #ifdef WSSMP_DBG
    1328                          if (ipp) chk(18, ipp, n);
    1329                          if (nnode) chk(19, nnode, n);
    1330 #endif
    1331                          if (ipp) snxt[ipp - 1] = nnode;
    1332                          else head[erscore[i - 1] - 1] = nnode;
    1333                          if (nnode) perm[nnode - 1] = ipp;
    1334                     }
    1335                }
    1336                currloc = 0;
    1337           }
    1338 #ifdef WSSMP_DBG
    1339           chk(20, node, n);
    1340 #endif
    1341           lsize[node - 1] = j5 - (xadj[node - 1] = j4);
    1342           dgree[node - 1] = nodeg;
    1343           if (fltag + n < 0 || fltag + n > MAXIW) {
    1344                for (i = 1; i <= n; i++) if (flag[i - 1]) flag[i - 1] = 1;
    1345                fltag = 2;
    1346           }
    1347           for (j3 = j4; j3 < j5; j3++) {
    1348 #ifdef WSSMP_DBG
    1349                chk(21, j3, adjln);
    1350 #endif
    1351                i = adjncy[j3 - 1];
    1352 #ifdef WSSMP_DBG
    1353                chk(22, i, n);
    1354 #endif
    1355                j = invp[i - 1];
    1356                if (j > 0) {
    1357                     i4 = fltag - (numii = -varbl[i - 1]);
    1358                     i2 = xadj[i - 1] + j;
    1359                     for (l = xadj[i - 1]; l < i2; l++) {
    1360 #ifdef WSSMP_DBG
    1361                          chk(23, l, adjln);
    1362 #endif
    1363                          tn = adjncy[l - 1];
    1364 #ifdef WSSMP_DBG
    1365                          chk(24, tn, n);
    1366 #endif
    1367                          j9 = flag[tn - 1];
    1368                          if (j9 >= fltag) j9 -= numii;
    1369                          else if (j9) j9 = dgree[tn - 1] + i4;
    1370                          flag[tn - 1] = j9;
    1371                     }
    1372                }
    1373           }
    1374           for (j3 = j4; j3 < j5; j3++) {
    1375 #ifdef WSSMP_DBG
    1376                chk(25, j3, adjln);
    1377 #endif
    1378                i = adjncy[j3 - 1];
    1379                i5 = deg = 0;
    1380 #ifdef WSSMP_DBG
    1381                chk(26, i, n);
    1382 #endif
    1383                j1 = (i4 = xadj[i - 1]) + invp[i - 1];
    1384                for (l = j0 = i4; l < j1; l++) {
    1385 #ifdef WSSMP_DBG
    1386                     chk(27, l, adjln);
    1387 #endif
    1388                     tn = adjncy[l - 1];
    1389 #ifdef WSSMP_DBG
    1390                     chk(70, tn, n);
    1391 #endif
    1392                     j8 = flag[tn - 1];
    1393                     if (j8) {
    1394                          deg += (j8 - fltag);
    1395 #ifdef WSSMP_DBG
    1396                          chk(28, i4, adjln);
    1397 #endif
    1398                          adjncy[i4-1] = tn;
    1399                          i5 += tn;
    1400                          i4++;
    1401                          while (i5 >= nm1) i5 -= nm1;
    1402                     }
    1403                }
    1404                invp[i - 1] = (j2 = i4) - j0 + 1;
    1405                i2 = j0 + lsize[i - 1];
    1406                for (l = j1; l < i2; l++) {
    1407 #ifdef WSSMP_DBG
    1408                     chk(29, l, adjln);
    1409 #endif
    1410                     j = adjncy[l - 1];
    1411 #ifdef WSSMP_DBG
    1412                     chk(30, j, n);
    1413 #endif
    1414                     numii = varbl[j - 1];
    1415                     if (numii > 0) {
    1416                          deg += numii;
    1417 #ifdef WSSMP_DBG
    1418                          chk(31, i4, adjln);
    1419 #endif
    1420                          adjncy[i4-1] = j;
    1421                          i5 += j;
    1422                          i4++;
    1423                          while (i5 >= nm1) i5 -= nm1;
    1424                     }
    1425                }
    1426                if (invp[i - 1] == 1 && j2 == i4) {
    1427                     numii = -varbl[i - 1];
    1428                     xadj[i - 1] = -node;
    1429                     nodeg -= numii;
    1430                     counter += numii;
    1431                     numpiv += numii;
    1432                     invp[i - 1] = varbl[i - 1] = 0;
    1433                } else {
    1434                     if (dgree[i - 1] > deg) dgree[i - 1] = deg;
    1435 #ifdef WSSMP_DBG
    1436                     chk(32, j2, adjln);
    1437                     chk(33, j0, adjln);
    1438 #endif
    1439                     adjncy[i4 - 1] = adjncy[j2 - 1];
    1440                     adjncy[j2 - 1] = adjncy[j0 - 1];
    1441                     adjncy[j0 - 1] = node;
    1442                     lsize[i - 1] = i4 - j0 + 1;
    1443                     i5++;
    1444 #ifdef WSSMP_DBG
    1445                     chk(35, i5, n);
    1446 #endif
    1447                     j = head[i5 - 1];
    1448                     if (j > 0) {
    1449                          snxt[i - 1] = perm[j - 1];
    1450                          perm[j - 1] = i;
    1451                     } else {
    1452                          snxt[i - 1] = -j;
    1453                          head[i5 - 1] = -i;
    1454                     }
    1455                     perm[i - 1] = i5;
    1456                }
    1457           }
    1458 #ifdef WSSMP_DBG
    1459           chk(36, node, n);
    1460 #endif
    1461           dgree[node - 1] = nodeg;
    1462           if (maxmum < nodeg) maxmum = nodeg;
    1463           fltag += maxmum;
    1464 #ifdef WSSMP_DBG
    1465           if (fltag + n + 1 < 0) printf("ERR**: fltag = %d (A)\n", fltag);
    1466 #endif
    1467           for (j3 = j4; j3 < j5; ++j3) {
    1468 #ifdef WSSMP_DBG
    1469                chk(37, j3, adjln);
    1470 #endif
    1471                i = adjncy[j3 - 1];
    1472 #ifdef WSSMP_DBG
    1473                chk(38, i, n);
    1474 #endif
    1475                if (varbl[i - 1] < 0) {
    1476                     i5 = perm[i - 1];
    1477 #ifdef WSSMP_DBG
    1478                     chk(39, i5, n);
    1479 #endif
    1480                     j = head[i5 - 1];
    1481                     if (j) {
    1482                          if (j < 0) {
    1483                               head[i5 - 1] = 0, i = -j;
    1484                          } else {
    1485 #ifdef WSSMP_DBG
    1486                               chk(40, j, n);
    1487 #endif
    1488                               i = perm[j - 1];
    1489                               perm[j - 1] = 0;
    1490                          }
    1491                          while (i) {
    1492 #ifdef WSSMP_DBG
    1493                               chk(41, i, n);
    1494 #endif
    1495                               if (!snxt[i - 1]) {
    1496                                    i = 0;
    1497                               } else {
    1498                                    k = invp[i - 1];
    1499                                    i2 = xadj[i - 1] + (scln = lsize[i - 1]);
    1500                                    for (l = xadj[i - 1] + 1; l < i2; l++) {
    1501 #ifdef WSSMP_DBG
    1502                                         chk(42, l, adjln);
    1503                                         chk(43, adjncy[l - 1], n);
    1504 #endif
    1505                                         flag[adjncy[l - 1] - 1] = fltag;
    1506                                    }
    1507                                    jpp = i;
    1508                                    j = snxt[i - 1];
    1509                                    while (j) {
    1510 #ifdef WSSMP_DBG
    1511                                         chk(44, j, n);
    1512 #endif
    1513                                         if (lsize[j - 1] == scln && invp[j - 1] == k) {
    1514                                              i2 = xadj[j - 1] + scln;
    1515                                              j8 = 1;
    1516                                              for (l = xadj[j - 1] + 1; l < i2; l++) {
    1517 #ifdef WSSMP_DBG
    1518                                                   chk(45, l, adjln);
    1519                                                   chk(46, adjncy[l - 1], n);
    1520 #endif
    1521                                                   if (flag[adjncy[l - 1] - 1] != fltag) {
    1522                                                        j8 = 0;
    1523                                                        break;
    1524                                                   }
    1525                                              }
    1526                                              if (j8) {
    1527                                                   xadj[j - 1] = -i;
    1528                                                   varbl[i - 1] += varbl[j - 1];
    1529                                                   varbl[j - 1] = invp[j - 1] = 0;
    1530 #ifdef WSSMP_DBG
    1531                                                   chk(47, j, n);
    1532                                                   chk(48, jpp, n);
    1533 #endif
    1534                                                   j = snxt[j - 1];
    1535                                                   snxt[jpp - 1] = j;
    1536                                              } else {
    1537                                                   jpp = j;
    1538 #ifdef WSSMP_DBG
    1539                                                   chk(49, j, n);
    1540 #endif
    1541                                                   j = snxt[j - 1];
    1542                                              }
    1543                                         } else {
    1544                                              jpp = j;
    1545 #ifdef WSSMP_DBG
    1546                                              chk(50, j, n);
    1547 #endif
    1548                                              j = snxt[j - 1];
    1549                                         }
    1550                                    }
    1551                                    fltag++;
    1552 #ifdef WSSMP_DBG
    1553                                    if (fltag + n + 1 < 0) printf("ERR**: fltag = %d (B)\n", fltag);
    1554 #endif
    1555 #ifdef WSSMP_DBG
    1556                                    chk(51, i, n);
    1557 #endif
    1558                                    i = snxt[i - 1];
    1559                               }
    1560                          }
    1561                     }
    1562                }
    1563           }
    1564           j8 = nm1 - counter;
    1565           switch (speed) {
    1566           case 1:
    1567                for (j = j3 = j4; j3 < j5; j3++) {
    1568 #ifdef WSSMP_DBG
    1569                     chk(52, j3, adjln);
    1570 #endif
    1571                     i = adjncy[j3 - 1];
    1572 #ifdef WSSMP_DBG
    1573                     chk(53, i, n);
    1574 #endif
    1575                     numii = varbl[i - 1];
    1576                     if (numii < 0) {
    1577                          k = 0;
    1578                          i4 = (l = xadj[i - 1]) + invp[i - 1];
    1579                          for (; l < i4; l++) {
    1580 #ifdef WSSMP_DBG
    1581                               chk(54, l, adjln);
    1582                               chk(55, adjncy[l - 1], n);
    1583 #endif
    1584                               i5 = dgree[adjncy[l - 1] - 1];
    1585                               if (k < i5) k = i5;
    1586                          }
    1587                          x = static_cast<double> (k - 1);
    1588                          varbl[i - 1] = abs(numii);
    1589                          j9 = dgree[i - 1] + nodeg;
    1590                          deg = (j8 > j9 ? j9 : j8) + numii;
    1591                          if (deg < 1) deg = 1;
    1592                          y = static_cast<double> (deg);
    1593                          dgree[i - 1] = deg;
    1594                          /*
     1549                  }
     1550                  if (j8) {
     1551                    xadj[j - 1] = -i;
     1552                    varbl[i - 1] += varbl[j - 1];
     1553                    varbl[j - 1] = invp[j - 1] = 0;
     1554#ifdef WSSMP_DBG
     1555                    chk(47, j, n);
     1556                    chk(48, jpp, n);
     1557#endif
     1558                    j = snxt[j - 1];
     1559                    snxt[jpp - 1] = j;
     1560                  } else {
     1561                    jpp = j;
     1562#ifdef WSSMP_DBG
     1563                    chk(49, j, n);
     1564#endif
     1565                    j = snxt[j - 1];
     1566                  }
     1567                } else {
     1568                  jpp = j;
     1569#ifdef WSSMP_DBG
     1570                  chk(50, j, n);
     1571#endif
     1572                  j = snxt[j - 1];
     1573                }
     1574              }
     1575              fltag++;
     1576#ifdef WSSMP_DBG
     1577              if (fltag + n + 1 < 0)
     1578                printf("ERR**: fltag = %d (B)\n", fltag);
     1579#endif
     1580#ifdef WSSMP_DBG
     1581              chk(51, i, n);
     1582#endif
     1583              i = snxt[i - 1];
     1584            }
     1585          }
     1586        }
     1587      }
     1588    }
     1589    j8 = nm1 - counter;
     1590    switch (speed) {
     1591    case 1:
     1592      for (j = j3 = j4; j3 < j5; j3++) {
     1593#ifdef WSSMP_DBG
     1594        chk(52, j3, adjln);
     1595#endif
     1596        i = adjncy[j3 - 1];
     1597#ifdef WSSMP_DBG
     1598        chk(53, i, n);
     1599#endif
     1600        numii = varbl[i - 1];
     1601        if (numii < 0) {
     1602          k = 0;
     1603          i4 = (l = xadj[i - 1]) + invp[i - 1];
     1604          for (; l < i4; l++) {
     1605#ifdef WSSMP_DBG
     1606            chk(54, l, adjln);
     1607            chk(55, adjncy[l - 1], n);
     1608#endif
     1609            i5 = dgree[adjncy[l - 1] - 1];
     1610            if (k < i5)
     1611              k = i5;
     1612          }
     1613          x = static_cast< double >(k - 1);
     1614          varbl[i - 1] = abs(numii);
     1615          j9 = dgree[i - 1] + nodeg;
     1616          deg = (j8 > j9 ? j9 : j8) + numii;
     1617          if (deg < 1)
     1618            deg = 1;
     1619          y = static_cast< double >(deg);
     1620          dgree[i - 1] = deg;
     1621          /*
    15951622                                    printf("%i %f should >= %i %f\n",deg,y,k-1,x);
    15961623                                    if (y < x) printf("** \n"); else printf("\n");
    15971624                         */
    1598                          y = y * y - y;
    1599                          x = y - (x * x - x);
    1600                          if (x < 1.1) x = 1.1;
    1601                          deg = static_cast<WSI>(sqrt(x));
    1602                          /*         if (deg < 1) deg = 1; */
    1603                          erscore[i - 1] = deg;
    1604 #ifdef WSSMP_DBG
    1605                          chk(56, deg, n - 1);
    1606 #endif
    1607                          nnode = head[deg - 1];
    1608                          if (nnode) perm[nnode - 1] = i;
    1609                          snxt[i - 1] = nnode;
    1610                          perm[i - 1] = 0;
    1611 #ifdef WSSMP_DBG
    1612                          chk(57, j, adjln);
    1613 #endif
    1614                          head[deg - 1] = adjncy[j-1] = i;
    1615                          j++;
    1616                          if (deg < mindeg) mindeg = deg;
    1617                     }
    1618                }
    1619                break;
    1620           case 2:
    1621                for (j = j3 = j4; j3 < j5; j3++) {
    1622 #ifdef WSSMP_DBG
    1623                     chk(58, j3, adjln);
    1624 #endif
    1625                     i = adjncy[j3 - 1];
    1626 #ifdef WSSMP_DBG
    1627                     chk(59, i, n);
    1628 #endif
    1629                     numii = varbl[i - 1];
    1630                     if (numii < 0) {
    1631                          x = static_cast<double> (dgree[adjncy[xadj[i - 1] - 1] - 1] - 1);
    1632                          varbl[i - 1] = abs(numii);
    1633                          j9 = dgree[i - 1] + nodeg;
    1634                          deg = (j8 > j9 ? j9 : j8) + numii;
    1635                          if (deg < 1) deg = 1;
    1636                          y = static_cast<double> (deg);
    1637                          dgree[i - 1] = deg;
    1638                          /*
     1625          y = y * y - y;
     1626          x = y - (x * x - x);
     1627          if (x < 1.1)
     1628            x = 1.1;
     1629          deg = static_cast< WSI >(sqrt(x));
     1630          /*         if (deg < 1) deg = 1; */
     1631          erscore[i - 1] = deg;
     1632#ifdef WSSMP_DBG
     1633          chk(56, deg, n - 1);
     1634#endif
     1635          nnode = head[deg - 1];
     1636          if (nnode)
     1637            perm[nnode - 1] = i;
     1638          snxt[i - 1] = nnode;
     1639          perm[i - 1] = 0;
     1640#ifdef WSSMP_DBG
     1641          chk(57, j, adjln);
     1642#endif
     1643          head[deg - 1] = adjncy[j - 1] = i;
     1644          j++;
     1645          if (deg < mindeg)
     1646            mindeg = deg;
     1647        }
     1648      }
     1649      break;
     1650    case 2:
     1651      for (j = j3 = j4; j3 < j5; j3++) {
     1652#ifdef WSSMP_DBG
     1653        chk(58, j3, adjln);
     1654#endif
     1655        i = adjncy[j3 - 1];
     1656#ifdef WSSMP_DBG
     1657        chk(59, i, n);
     1658#endif
     1659        numii = varbl[i - 1];
     1660        if (numii < 0) {
     1661          x = static_cast< double >(dgree[adjncy[xadj[i - 1] - 1] - 1] - 1);
     1662          varbl[i - 1] = abs(numii);
     1663          j9 = dgree[i - 1] + nodeg;
     1664          deg = (j8 > j9 ? j9 : j8) + numii;
     1665          if (deg < 1)
     1666            deg = 1;
     1667          y = static_cast< double >(deg);
     1668          dgree[i - 1] = deg;
     1669          /*
    16391670                                    printf("%i %f should >= %i %f",deg,y,dgree[adjncy[xadj[i - 1] - 1] - 1]-1,x);
    16401671                                    if (y < x) printf("** \n"); else printf("\n");
    16411672                         */
    1642                          y = y * y - y;
    1643                          x = y - (x * x - x);
    1644                          if (x < 1.1) x = 1.1;
    1645                          deg = static_cast<WSI>(sqrt(x));
    1646                          /*         if (deg < 1) deg = 1; */
    1647                          erscore[i - 1] = deg;
    1648 #ifdef WSSMP_DBG
    1649                          chk(60, deg, n - 1);
    1650 #endif
    1651                          nnode = head[deg - 1];
    1652                          if (nnode) perm[nnode - 1] = i;
    1653                          snxt[i - 1] = nnode;
    1654                          perm[i - 1] = 0;
    1655 #ifdef WSSMP_DBG
    1656                          chk(61, j, adjln);
    1657 #endif
    1658                          head[deg - 1] = adjncy[j-1] = i;
    1659                          j++;
    1660                          if (deg < mindeg) mindeg = deg;
    1661                     }
    1662                }
    1663                break;
    1664           default:
    1665                for (j = j3 = j4; j3 < j5; j3++) {
    1666 #ifdef WSSMP_DBG
    1667                     chk(62, j3, adjln);
    1668 #endif
    1669                     i = adjncy[j3 - 1];
    1670 #ifdef WSSMP_DBG
    1671                     chk(63, i, n);
    1672 #endif
    1673                     numii = varbl[i - 1];
    1674                     if (numii < 0) {
    1675                          varbl[i - 1] = abs(numii);
    1676                          j9 = dgree[i - 1] + nodeg;
    1677                          deg = (j8 > j9 ? j9 : j8) + numii;
    1678                          if (deg < 1) deg = 1;
    1679                          dgree[i - 1] = deg;
    1680 #ifdef WSSMP_DBG
    1681                          chk(64, deg, n - 1);
    1682 #endif
    1683                          nnode = head[deg - 1];
    1684                          if (nnode) perm[nnode - 1] = i;
    1685                          snxt[i - 1] = nnode;
    1686                          perm[i - 1] = 0;
    1687 #ifdef WSSMP_DBG
    1688                          chk(65, j, adjln);
    1689 #endif
    1690                          head[deg - 1] = adjncy[j-1] = i;
    1691                          j++;
    1692                          if (deg < mindeg) mindeg = deg;
    1693                     }
    1694                }
    1695                break;
    1696           }
    1697           if (currloc) {
    1698 #ifdef WSSMP_DBG
    1699                chk(66, node, n);
    1700 #endif
    1701                locatns += (lsize[node - 1] - currloc), locaux = j;
    1702           }
    1703           varbl[node - 1] = numpiv + nodeg;
    1704           lsize[node - 1] = j - j4;
    1705           if (!lsize[node - 1]) flag[node - 1] = xadj[node - 1] = 0;
    1706      }
    1707      for (l = 1; l <= n; l++) {
    1708           if (!invp[l - 1]) {
    1709                for (i = -xadj[l - 1]; invp[i - 1] >= 0; i = -xadj[i - 1]) {};
    1710                tn = i;
    1711 #ifdef WSSMP_DBG
    1712                chk(67, tn, n);
    1713 #endif
    1714                k = -invp[tn - 1];
    1715                i = l;
    1716 #ifdef WSSMP_DBG
    1717                chk(68, i, n);
    1718 #endif
    1719                while (invp[i - 1] >= 0) {
    1720                     nnode = -xadj[i - 1];
    1721                     xadj[i - 1] = -tn;
    1722                     if (!invp[i - 1]) invp[i - 1] = k++;
    1723                     i = nnode;
    1724                }
    1725                invp[tn - 1] = -k;
    1726           }
    1727      }
    1728      for (l = 0; l < n; l++) {
    1729           i = abs(invp[l]);
    1730 #ifdef WSSMP_DBG
    1731           chk(69, i, n);
    1732 #endif
    1733           invp[l] = i;
    1734           perm[i - 1] = l + 1;
    1735      }
    1736      return;
     1673          y = y * y - y;
     1674          x = y - (x * x - x);
     1675          if (x < 1.1)
     1676            x = 1.1;
     1677          deg = static_cast< WSI >(sqrt(x));
     1678          /*         if (deg < 1) deg = 1; */
     1679          erscore[i - 1] = deg;
     1680#ifdef WSSMP_DBG
     1681          chk(60, deg, n - 1);
     1682#endif
     1683          nnode = head[deg - 1];
     1684          if (nnode)
     1685            perm[nnode - 1] = i;
     1686          snxt[i - 1] = nnode;
     1687          perm[i - 1] = 0;
     1688#ifdef WSSMP_DBG
     1689          chk(61, j, adjln);
     1690#endif
     1691          head[deg - 1] = adjncy[j - 1] = i;
     1692          j++;
     1693          if (deg < mindeg)
     1694            mindeg = deg;
     1695        }
     1696      }
     1697      break;
     1698    default:
     1699      for (j = j3 = j4; j3 < j5; j3++) {
     1700#ifdef WSSMP_DBG
     1701        chk(62, j3, adjln);
     1702#endif
     1703        i = adjncy[j3 - 1];
     1704#ifdef WSSMP_DBG
     1705        chk(63, i, n);
     1706#endif
     1707        numii = varbl[i - 1];
     1708        if (numii < 0) {
     1709          varbl[i - 1] = abs(numii);
     1710          j9 = dgree[i - 1] + nodeg;
     1711          deg = (j8 > j9 ? j9 : j8) + numii;
     1712          if (deg < 1)
     1713            deg = 1;
     1714          dgree[i - 1] = deg;
     1715#ifdef WSSMP_DBG
     1716          chk(64, deg, n - 1);
     1717#endif
     1718          nnode = head[deg - 1];
     1719          if (nnode)
     1720            perm[nnode - 1] = i;
     1721          snxt[i - 1] = nnode;
     1722          perm[i - 1] = 0;
     1723#ifdef WSSMP_DBG
     1724          chk(65, j, adjln);
     1725#endif
     1726          head[deg - 1] = adjncy[j - 1] = i;
     1727          j++;
     1728          if (deg < mindeg)
     1729            mindeg = deg;
     1730        }
     1731      }
     1732      break;
     1733    }
     1734    if (currloc) {
     1735#ifdef WSSMP_DBG
     1736      chk(66, node, n);
     1737#endif
     1738      locatns += (lsize[node - 1] - currloc), locaux = j;
     1739    }
     1740    varbl[node - 1] = numpiv + nodeg;
     1741    lsize[node - 1] = j - j4;
     1742    if (!lsize[node - 1])
     1743      flag[node - 1] = xadj[node - 1] = 0;
     1744  }
     1745  for (l = 1; l <= n; l++) {
     1746    if (!invp[l - 1]) {
     1747      for (i = -xadj[l - 1]; invp[i - 1] >= 0; i = -xadj[i - 1]) {
     1748      };
     1749      tn = i;
     1750#ifdef WSSMP_DBG
     1751      chk(67, tn, n);
     1752#endif
     1753      k = -invp[tn - 1];
     1754      i = l;
     1755#ifdef WSSMP_DBG
     1756      chk(68, i, n);
     1757#endif
     1758      while (invp[i - 1] >= 0) {
     1759        nnode = -xadj[i - 1];
     1760        xadj[i - 1] = -tn;
     1761        if (!invp[i - 1])
     1762          invp[i - 1] = k++;
     1763        i = nnode;
     1764      }
     1765      invp[tn - 1] = -k;
     1766    }
     1767  }
     1768  for (l = 0; l < n; l++) {
     1769    i = abs(invp[l]);
     1770#ifdef WSSMP_DBG
     1771    chk(69, i, n);
     1772#endif
     1773    invp[l] = i;
     1774    perm[i - 1] = l + 1;
     1775  }
     1776  return;
    17371777}
    17381778// This code is not needed, but left in in case needed sometime
     
    18191859#endif
    18201860// Orders rows
    1821 int
    1822 ClpCholeskyBase::orderAMD()
     1861int ClpCholeskyBase::orderAMD()
    18231862{
    1824      permuteInverse_ = new int [numberRows_];
    1825      permute_ = new int[numberRows_];
    1826      // Do ordering
    1827      int returnCode = 0;
    1828      // get full matrix
    1829      int space = 2 * sizeFactor_ + 10000 + 4 * numberRows_;
    1830      int * temp = new int [space];
    1831      int * count = new int [numberRows_];
    1832      int * tempStart = new int [numberRows_+1];
    1833      memset(count, 0, numberRows_ * sizeof(int));
    1834      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1835        count[iRow] += static_cast<int>(choleskyStart_[iRow+1] - choleskyStart_[iRow] - 1);
    1836           for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
    1837                int jRow = choleskyRow_[j];
    1838                count[jRow]++;
    1839           }
    1840      }
     1863  permuteInverse_ = new int[numberRows_];
     1864  permute_ = new int[numberRows_];
     1865  // Do ordering
     1866  int returnCode = 0;
     1867  // get full matrix
     1868  int space = 2 * sizeFactor_ + 10000 + 4 * numberRows_;
     1869  int *temp = new int[space];
     1870  int *count = new int[numberRows_];
     1871  int *tempStart = new int[numberRows_ + 1];
     1872  memset(count, 0, numberRows_ * sizeof(int));
     1873  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1874    count[iRow] += static_cast< int >(choleskyStart_[iRow + 1] - choleskyStart_[iRow] - 1);
     1875    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
     1876      int jRow = choleskyRow_[j];
     1877      count[jRow]++;
     1878    }
     1879  }
    18411880#define OFFSET 1
    1842      int sizeFactor = 0;
    1843      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1844           int length = count[iRow];
    1845           permute_[iRow] = length;
    1846           // add 1 to starts
    1847           tempStart[iRow] = sizeFactor + OFFSET;
    1848           count[iRow] = sizeFactor;
    1849           sizeFactor += length;
    1850      }
    1851      tempStart[numberRows_] = sizeFactor + OFFSET;
    1852      // add 1 to rows
    1853      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1854           assert (choleskyRow_[choleskyStart_[iRow]] == iRow);
    1855           for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
    1856                int jRow = choleskyRow_[j];
    1857                int put = count[iRow];
    1858                temp[put] = jRow + OFFSET;
    1859                count[iRow]++;
    1860                put = count[jRow];
    1861                temp[put] = iRow + OFFSET;
    1862                count[jRow]++;
    1863           }
    1864      }
    1865      for (int iRow = 1; iRow < numberRows_; iRow++)
    1866           assert (count[iRow-1] == tempStart[iRow] - OFFSET);
    1867      delete [] choleskyRow_;
    1868      choleskyRow_ = temp;
    1869      delete [] choleskyStart_;
    1870      choleskyStart_ = tempStart;
    1871      int locaux = sizeFactor + OFFSET;
    1872      delete [] count;
    1873      int speed = integerParameters_[0];
    1874      if (speed < 1 || speed > 2)
    1875           speed = 3;
    1876      int * use = new int [((speed<3) ? 7 : 6)*numberRows_];
    1877      int * dgree = use;
    1878      int * varbl = dgree + numberRows_;
    1879      int * snxt = varbl + numberRows_;
    1880      int * head = snxt + numberRows_;
    1881      int * lsize = head + numberRows_;
    1882      int * flag = lsize + numberRows_;
    1883      int * erscore;
    1884      for (int i = 0; i < numberRows_; i++) {
    1885           dgree[i] = choleskyStart_[i+1] - choleskyStart_[i];
    1886           head[i] = dgree[i];
    1887           snxt[i] = 0;
    1888           permute_[i] = 0;
    1889           permuteInverse_[i] = 0;
    1890           head[i] = 0;
    1891           flag[i] = 1;
    1892           varbl[i] = 1;
    1893           lsize[i] = dgree[i];
    1894      }
    1895      if (speed < 3) {
    1896           erscore = flag + numberRows_;
    1897           for (int i = 0; i < numberRows_; i++)
    1898                erscore[i] = dgree[i];
    1899      } else {
    1900           erscore = dgree;
    1901      }
    1902      myamlf(numberRows_, choleskyStart_, choleskyRow_,
    1903             dgree, varbl, snxt, permute_, permuteInverse_,
    1904             head, lsize, flag, erscore, locaux, space, speed);
    1905      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1906           permute_[iRow]--;
    1907      }
    1908      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1909           permuteInverse_[permute_[iRow]] = iRow;
    1910      }
    1911      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1912           assert (permuteInverse_[iRow] >= 0 && permuteInverse_[iRow] < numberRows_);
    1913      }
    1914      delete [] use;
    1915      delete [] choleskyRow_;
    1916      choleskyRow_ = NULL;
    1917      delete [] choleskyStart_;
    1918      choleskyStart_ = NULL;
    1919      return returnCode;
     1881  int sizeFactor = 0;
     1882  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1883    int length = count[iRow];
     1884    permute_[iRow] = length;
     1885    // add 1 to starts
     1886    tempStart[iRow] = sizeFactor + OFFSET;
     1887    count[iRow] = sizeFactor;
     1888    sizeFactor += length;
     1889  }
     1890  tempStart[numberRows_] = sizeFactor + OFFSET;
     1891  // add 1 to rows
     1892  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1893    assert(choleskyRow_[choleskyStart_[iRow]] == iRow);
     1894    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
     1895      int jRow = choleskyRow_[j];
     1896      int put = count[iRow];
     1897      temp[put] = jRow + OFFSET;
     1898      count[iRow]++;
     1899      put = count[jRow];
     1900      temp[put] = iRow + OFFSET;
     1901      count[jRow]++;
     1902    }
     1903  }
     1904  for (int iRow = 1; iRow < numberRows_; iRow++)
     1905    assert(count[iRow - 1] == tempStart[iRow] - OFFSET);
     1906  delete[] choleskyRow_;
     1907  choleskyRow_ = temp;
     1908  delete[] choleskyStart_;
     1909  choleskyStart_ = tempStart;
     1910  int locaux = sizeFactor + OFFSET;
     1911  delete[] count;
     1912  int speed = integerParameters_[0];
     1913  if (speed < 1 || speed > 2)
     1914    speed = 3;
     1915  int *use = new int[((speed < 3) ? 7 : 6) * numberRows_];
     1916  int *dgree = use;
     1917  int *varbl = dgree + numberRows_;
     1918  int *snxt = varbl + numberRows_;
     1919  int *head = snxt + numberRows_;
     1920  int *lsize = head + numberRows_;
     1921  int *flag = lsize + numberRows_;
     1922  int *erscore;
     1923  for (int i = 0; i < numberRows_; i++) {
     1924    dgree[i] = choleskyStart_[i + 1] - choleskyStart_[i];
     1925    head[i] = dgree[i];
     1926    snxt[i] = 0;
     1927    permute_[i] = 0;
     1928    permuteInverse_[i] = 0;
     1929    head[i] = 0;
     1930    flag[i] = 1;
     1931    varbl[i] = 1;
     1932    lsize[i] = dgree[i];
     1933  }
     1934  if (speed < 3) {
     1935    erscore = flag + numberRows_;
     1936    for (int i = 0; i < numberRows_; i++)
     1937      erscore[i] = dgree[i];
     1938  } else {
     1939    erscore = dgree;
     1940  }
     1941  myamlf(numberRows_, choleskyStart_, choleskyRow_,
     1942    dgree, varbl, snxt, permute_, permuteInverse_,
     1943    head, lsize, flag, erscore, locaux, space, speed);
     1944  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1945    permute_[iRow]--;
     1946  }
     1947  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1948    permuteInverse_[permute_[iRow]] = iRow;
     1949  }
     1950  for (int iRow = 0; iRow < numberRows_; iRow++) {
     1951    assert(permuteInverse_[iRow] >= 0 && permuteInverse_[iRow] < numberRows_);
     1952  }
     1953  delete[] use;
     1954  delete[] choleskyRow_;
     1955  choleskyRow_ = NULL;
     1956  delete[] choleskyStart_;
     1957  choleskyStart_ = NULL;
     1958  return returnCode;
    19201959}
    19211960/* Does Symbolic factorization given permutation.
     
    19231962   user must provide factorize and solve.  Otherwise the default factorization is used
    19241963   returns non-zero if not enough memory */
    1925 int
    1926 ClpCholeskyBase::symbolic()
     1964int ClpCholeskyBase::symbolic()
    19271965{
    1928      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    1929      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    1930      const int * row = model_->clpMatrix()->getIndices();
    1931      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    1932      const int * rowLength = rowCopy_->getVectorLengths();
    1933      const int * column = rowCopy_->getIndices();
    1934      int numberRowsModel = model_->numberRows();
    1935      int numberColumns = model_->numberColumns();
    1936      int numberTotal = numberColumns + numberRowsModel;
    1937      CoinPackedMatrix * quadratic = NULL;
    1938      ClpQuadraticObjective * quadraticObj =
    1939           (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
    1940      if (quadraticObj)
    1941           quadratic = quadraticObj->quadraticObjective();
    1942      // We need an array for counts
    1943      int * used = new int[numberRows_+1];
    1944      // If KKT then re-order so negative first
    1945      if (doKKT_) {
    1946           int nn = 0;
    1947           int np = 0;
    1948           int iRow;
     1966  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     1967  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     1968  const int *row = model_->clpMatrix()->getIndices();
     1969  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     1970  const int *rowLength = rowCopy_->getVectorLengths();
     1971  const int *column = rowCopy_->getIndices();
     1972  int numberRowsModel = model_->numberRows();
     1973  int numberColumns = model_->numberColumns();
     1974  int numberTotal = numberColumns + numberRowsModel;
     1975  CoinPackedMatrix *quadratic = NULL;
     1976  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
     1977  if (quadraticObj)
     1978    quadratic = quadraticObj->quadraticObjective();
     1979  // We need an array for counts
     1980  int *used = new int[numberRows_ + 1];
     1981  // If KKT then re-order so negative first
     1982  if (doKKT_) {
     1983    int nn = 0;
     1984    int np = 0;
     1985    int iRow;
     1986    for (iRow = 0; iRow < numberRows_; iRow++) {
     1987      int originalRow = permute_[iRow];
     1988      if (originalRow < numberTotal)
     1989        permute_[nn++] = originalRow;
     1990      else
     1991        used[np++] = originalRow;
     1992    }
     1993    CoinMemcpyN(used, np, permute_ + nn);
     1994    for (iRow = 0; iRow < numberRows_; iRow++)
     1995      permuteInverse_[permute_[iRow]] = iRow;
     1996  }
     1997  CoinZeroN(used, numberRows_);
     1998  int iRow;
     1999  int iColumn;
     2000  bool noMemory = false;
     2001  int *Astart = new int[numberRows_ + 1];
     2002  int *Arow = NULL;
     2003  try {
     2004    Arow = new int[sizeFactor_];
     2005  } catch (...) {
     2006    // no memory
     2007    delete[] Astart;
     2008    return -1;
     2009  }
     2010  choleskyStart_ = new int[numberRows_ + 1];
     2011  link_ = new int[numberRows_];
     2012  workInteger_ = new int[numberRows_];
     2013  indexStart_ = new int[numberRows_];
     2014  clique_ = new int[numberRows_];
     2015  // Redo so permuted upper triangular
     2016  sizeFactor_ = 0;
     2017  int *which = Arow;
     2018  if (!doKKT_) {
     2019    for (iRow = 0; iRow < numberRows_; iRow++) {
     2020      int number = 0;
     2021      int iOriginalRow = permute_[iRow];
     2022      Astart[iRow] = sizeFactor_;
     2023      CoinBigIndex startRow = rowStart[iOriginalRow];
     2024      CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
     2025      for (CoinBigIndex k = startRow; k < endRow; k++) {
     2026        int iColumn = column[k];
     2027        if (!whichDense_ || !whichDense_[iColumn]) {
     2028          CoinBigIndex start = columnStart[iColumn];
     2029          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2030          for (CoinBigIndex j = start; j < end; j++) {
     2031            int jRow = row[j];
     2032            int jNewRow = permuteInverse_[jRow];
     2033            if (jNewRow < iRow) {
     2034              if (!used[jNewRow]) {
     2035                used[jNewRow] = 1;
     2036                which[number++] = jNewRow;
     2037              }
     2038            }
     2039          }
     2040        }
     2041      }
     2042      sizeFactor_ += number;
     2043      int j;
     2044      for (j = 0; j < number; j++)
     2045        used[which[j]] = 0;
     2046      // Sort
     2047      std::sort(which, which + number);
     2048      // move which on
     2049      which += number;
     2050    }
     2051  } else {
     2052    // KKT
     2053    // transpose
     2054    ClpMatrixBase *rowCopy = model_->clpMatrix()->reverseOrderedCopy();
     2055    const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
     2056    const int *rowLength = rowCopy->getVectorLengths();
     2057    const int *column = rowCopy->getIndices();
     2058    // temp
     2059    bool permuted = false;
     2060    for (iRow = 0; iRow < numberRows_; iRow++) {
     2061      if (permute_[iRow] != iRow) {
     2062        permuted = true;
     2063        break;
     2064      }
     2065    }
     2066    if (permuted) {
     2067      // Need to permute - ugly
     2068      if (!quadratic) {
     2069        for (iRow = 0; iRow < numberRows_; iRow++) {
     2070          Astart[iRow] = sizeFactor_;
     2071          int iOriginalRow = permute_[iRow];
     2072          if (iOriginalRow < numberColumns) {
     2073            // A may be upper triangular by mistake
     2074            iColumn = iOriginalRow;
     2075            CoinBigIndex start = columnStart[iColumn];
     2076            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2077            for (CoinBigIndex j = start; j < end; j++) {
     2078              int kRow = row[j] + numberTotal;
     2079              kRow = permuteInverse_[kRow];
     2080              if (kRow < iRow)
     2081                Arow[sizeFactor_++] = kRow;
     2082            }
     2083          } else if (iOriginalRow < numberTotal) {
     2084            int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
     2085            if (kRow < iRow)
     2086              Arow[sizeFactor_++] = kRow;
     2087          } else {
     2088            int kRow = iOriginalRow - numberTotal;
     2089            CoinBigIndex start = rowStart[kRow];
     2090            CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
     2091            for (CoinBigIndex j = start; j < end; j++) {
     2092              int jRow = column[j];
     2093              int jNewRow = permuteInverse_[jRow];
     2094              if (jNewRow < iRow)
     2095                Arow[sizeFactor_++] = jNewRow;
     2096            }
     2097            // slack - should it be permute
     2098            kRow = permuteInverse_[kRow + numberColumns];
     2099            if (kRow < iRow)
     2100              Arow[sizeFactor_++] = kRow;
     2101          }
     2102          // Sort
     2103          std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
     2104        }
     2105      } else {
     2106        // quadratic
     2107        // transpose
     2108        CoinPackedMatrix quadraticT;
     2109        quadraticT.reverseOrderedCopyOf(*quadratic);
     2110        const int *columnQuadratic = quadraticT.getIndices();
     2111        const CoinBigIndex *columnQuadraticStart = quadraticT.getVectorStarts();
     2112        const int *columnQuadraticLength = quadraticT.getVectorLengths();
     2113        for (iRow = 0; iRow < numberRows_; iRow++) {
     2114          Astart[iRow] = sizeFactor_;
     2115          int iOriginalRow = permute_[iRow];
     2116          if (iOriginalRow < numberColumns) {
     2117            // Quadratic bit
     2118            CoinBigIndex j;
     2119            for (j = columnQuadraticStart[iOriginalRow];
     2120                 j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
     2121              int jColumn = columnQuadratic[j];
     2122              int jNewColumn = permuteInverse_[jColumn];
     2123              if (jNewColumn < iRow)
     2124                Arow[sizeFactor_++] = jNewColumn;
     2125            }
     2126            // A may be upper triangular by mistake
     2127            iColumn = iOriginalRow;
     2128            CoinBigIndex start = columnStart[iColumn];
     2129            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2130            for (j = start; j < end; j++) {
     2131              int kRow = row[j] + numberTotal;
     2132              kRow = permuteInverse_[kRow];
     2133              if (kRow < iRow)
     2134                Arow[sizeFactor_++] = kRow;
     2135            }
     2136          } else if (iOriginalRow < numberTotal) {
     2137            int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
     2138            if (kRow < iRow)
     2139              Arow[sizeFactor_++] = kRow;
     2140          } else {
     2141            int kRow = iOriginalRow - numberTotal;
     2142            CoinBigIndex start = rowStart[kRow];
     2143            CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
     2144            for (CoinBigIndex j = start; j < end; j++) {
     2145              int jRow = column[j];
     2146              int jNewRow = permuteInverse_[jRow];
     2147              if (jNewRow < iRow)
     2148                Arow[sizeFactor_++] = jNewRow;
     2149            }
     2150            // slack - should it be permute
     2151            kRow = permuteInverse_[kRow + numberColumns];
     2152            if (kRow < iRow)
     2153              Arow[sizeFactor_++] = kRow;
     2154          }
     2155          // Sort
     2156          std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
     2157        }
     2158      }
     2159    } else {
     2160      if (!quadratic) {
     2161        for (iRow = 0; iRow < numberRows_; iRow++) {
     2162          Astart[iRow] = sizeFactor_;
     2163        }
     2164      } else {
     2165        // Quadratic
     2166        // transpose
     2167        CoinPackedMatrix quadraticT;
     2168        quadraticT.reverseOrderedCopyOf(*quadratic);
     2169        const int *columnQuadratic = quadraticT.getIndices();
     2170        const CoinBigIndex *columnQuadraticStart = quadraticT.getVectorStarts();
     2171        const int *columnQuadraticLength = quadraticT.getVectorLengths();
     2172        //const double * quadraticElement = quadraticT.getElements();
     2173        for (iRow = 0; iRow < numberColumns; iRow++) {
     2174          int iOriginalRow = permute_[iRow];
     2175          Astart[iRow] = sizeFactor_;
     2176          for (CoinBigIndex j = columnQuadraticStart[iOriginalRow];
     2177               j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
     2178            int jColumn = columnQuadratic[j];
     2179            int jNewColumn = permuteInverse_[jColumn];
     2180            if (jNewColumn < iRow)
     2181              Arow[sizeFactor_++] = jNewColumn;
     2182          }
     2183        }
     2184      }
     2185      int iRow;
     2186      // slacks
     2187      for (iRow = 0; iRow < numberRowsModel; iRow++) {
     2188        Astart[iRow + numberColumns] = sizeFactor_;
     2189      }
     2190      for (iRow = 0; iRow < numberRowsModel; iRow++) {
     2191        Astart[iRow + numberTotal] = sizeFactor_;
     2192        CoinBigIndex start = rowStart[iRow];
     2193        CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
     2194        for (CoinBigIndex j = start; j < end; j++) {
     2195          Arow[sizeFactor_++] = column[j];
     2196        }
     2197        // slack
     2198        Arow[sizeFactor_++] = numberColumns + iRow;
     2199      }
     2200    }
     2201    delete rowCopy;
     2202  }
     2203  Astart[numberRows_] = sizeFactor_;
     2204  firstDense_ = numberRows_;
     2205  symbolic1(Astart, Arow);
     2206  // Now fill in indices
     2207  try {
     2208    // too big
     2209    choleskyRow_ = new int[sizeFactor_];
     2210  } catch (...) {
     2211    // no memory
     2212    noMemory = true;
     2213  }
     2214  double sizeFactor = sizeFactor_;
     2215  if (!noMemory) {
     2216    // Do lower triangular
     2217    sizeFactor_ = 0;
     2218    int *which = Arow;
     2219    if (!doKKT_) {
     2220      for (iRow = 0; iRow < numberRows_; iRow++) {
     2221        int number = 0;
     2222        int iOriginalRow = permute_[iRow];
     2223        Astart[iRow] = sizeFactor_;
     2224        if (!rowsDropped_[iOriginalRow]) {
     2225          CoinBigIndex startRow = rowStart[iOriginalRow];
     2226          CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
     2227          for (CoinBigIndex k = startRow; k < endRow; k++) {
     2228            int iColumn = column[k];
     2229            if (!whichDense_ || !whichDense_[iColumn]) {
     2230              CoinBigIndex start = columnStart[iColumn];
     2231              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2232              for (CoinBigIndex j = start; j < end; j++) {
     2233                int jRow = row[j];
     2234                int jNewRow = permuteInverse_[jRow];
     2235                if (jNewRow > iRow && !rowsDropped_[jRow]) {
     2236                  if (!used[jNewRow]) {
     2237                    used[jNewRow] = 1;
     2238                    which[number++] = jNewRow;
     2239                  }
     2240                }
     2241              }
     2242            }
     2243          }
     2244          sizeFactor_ += number;
     2245          int j;
     2246          for (j = 0; j < number; j++)
     2247            used[which[j]] = 0;
     2248          // Sort
     2249          std::sort(which, which + number);
     2250          // move which on
     2251          which += number;
     2252        }
     2253      }
     2254    } else {
     2255      // KKT
     2256      // temp
     2257      bool permuted = false;
     2258      for (iRow = 0; iRow < numberRows_; iRow++) {
     2259        if (permute_[iRow] != iRow) {
     2260          permuted = true;
     2261          break;
     2262        }
     2263      }
     2264      // but fake it
     2265      for (iRow = 0; iRow < numberRows_; iRow++) {
     2266        //permute_[iRow]=iRow; // force no permute
     2267        //permuteInverse_[permute_[iRow]]=iRow;
     2268      }
     2269      if (permuted) {
     2270        // Need to permute - ugly
     2271        if (!quadratic) {
    19492272          for (iRow = 0; iRow < numberRows_; iRow++) {
    1950                int originalRow = permute_[iRow];
    1951                if (originalRow < numberTotal)
    1952                     permute_[nn++] = originalRow;
    1953                else
    1954                     used[np++] = originalRow;
    1955           }
    1956           CoinMemcpyN(used, np, permute_ + nn);
    1957           for (iRow = 0; iRow < numberRows_; iRow++)
    1958                permuteInverse_[permute_[iRow]] = iRow;
    1959      }
    1960      CoinZeroN(used, numberRows_);
    1961      int iRow;
    1962      int iColumn;
    1963      bool noMemory = false;
    1964      int * Astart = new int[numberRows_+1];
    1965      int * Arow = NULL;
    1966      try {
    1967           Arow = new int [sizeFactor_];
    1968      } catch (...) {
    1969           // no memory
    1970           delete [] Astart;
    1971           return -1;
    1972      }
    1973      choleskyStart_ = new int[numberRows_+1];
    1974      link_ = new int[numberRows_];
    1975      workInteger_ = new int [numberRows_];
    1976      indexStart_ = new int[numberRows_];
    1977      clique_ = new int[numberRows_];
    1978      // Redo so permuted upper triangular
    1979      sizeFactor_ = 0;
    1980      int * which = Arow;
    1981      if (!doKKT_) {
     2273            Astart[iRow] = sizeFactor_;
     2274            int iOriginalRow = permute_[iRow];
     2275            if (iOriginalRow < numberColumns) {
     2276              iColumn = iOriginalRow;
     2277              CoinBigIndex start = columnStart[iColumn];
     2278              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2279              for (CoinBigIndex j = start; j < end; j++) {
     2280                int kRow = row[j] + numberTotal;
     2281                kRow = permuteInverse_[kRow];
     2282                if (kRow > iRow)
     2283                  Arow[sizeFactor_++] = kRow;
     2284              }
     2285            } else if (iOriginalRow < numberTotal) {
     2286              int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
     2287              if (kRow > iRow)
     2288                Arow[sizeFactor_++] = kRow;
     2289            } else {
     2290              int kRow = iOriginalRow - numberTotal;
     2291              CoinBigIndex start = rowStart[kRow];
     2292              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
     2293              for (CoinBigIndex j = start; j < end; j++) {
     2294                int jRow = column[j];
     2295                int jNewRow = permuteInverse_[jRow];
     2296                if (jNewRow > iRow)
     2297                  Arow[sizeFactor_++] = jNewRow;
     2298              }
     2299              // slack - should it be permute
     2300              kRow = permuteInverse_[kRow + numberColumns];
     2301              if (kRow > iRow)
     2302                Arow[sizeFactor_++] = kRow;
     2303            }
     2304            // Sort
     2305            std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
     2306          }
     2307        } else {
     2308          // quadratic
     2309          const int *columnQuadratic = quadratic->getIndices();
     2310          const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
     2311          const int *columnQuadraticLength = quadratic->getVectorLengths();
    19822312          for (iRow = 0; iRow < numberRows_; iRow++) {
    1983                int number = 0;
    1984                int iOriginalRow = permute_[iRow];
    1985                Astart[iRow] = sizeFactor_;
    1986                CoinBigIndex startRow = rowStart[iOriginalRow];
    1987                CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
    1988                for (CoinBigIndex k = startRow; k < endRow; k++) {
    1989                     int iColumn = column[k];
    1990                     if (!whichDense_ || !whichDense_[iColumn]) {
    1991                          CoinBigIndex start = columnStart[iColumn];
    1992                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    1993                          for (CoinBigIndex j = start; j < end; j++) {
    1994                               int jRow = row[j];
    1995                               int jNewRow = permuteInverse_[jRow];
    1996                               if (jNewRow < iRow) {
    1997                                    if (!used[jNewRow]) {
    1998                                         used[jNewRow] = 1;
    1999                                         which[number++] = jNewRow;
    2000                                    }
    2001                               }
    2002                          }
    2003                     }
    2004                }
    2005                sizeFactor_ += number;
    2006                int j;
    2007                for (j = 0; j < number; j++)
    2008                     used[which[j]] = 0;
    2009                // Sort
    2010                std::sort(which, which + number);
    2011                // move which on
    2012                which += number;
    2013           }
    2014      } else {
    2015           // KKT
    2016           // transpose
    2017           ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
    2018           const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    2019           const int * rowLength = rowCopy->getVectorLengths();
    2020           const int * column = rowCopy->getIndices();
    2021           // temp
    2022           bool permuted = false;
    2023           for (iRow = 0; iRow < numberRows_; iRow++) {
    2024                if (permute_[iRow] != iRow) {
    2025                     permuted = true;
    2026                     break;
    2027                }
    2028           }
    2029           if (permuted) {
    2030                // Need to permute - ugly
    2031                if (!quadratic) {
    2032                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2033                          Astart[iRow] = sizeFactor_;
    2034                          int iOriginalRow = permute_[iRow];
    2035                          if (iOriginalRow < numberColumns) {
    2036                               // A may be upper triangular by mistake
    2037                               iColumn = iOriginalRow;
    2038                               CoinBigIndex start = columnStart[iColumn];
    2039                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2040                               for (CoinBigIndex j = start; j < end; j++) {
    2041                                    int kRow = row[j] + numberTotal;
    2042                                    kRow = permuteInverse_[kRow];
    2043                                    if (kRow < iRow)
    2044                                         Arow[sizeFactor_++] = kRow;
    2045                               }
    2046                          } else if (iOriginalRow < numberTotal) {
    2047                               int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
    2048                               if (kRow < iRow)
    2049                                    Arow[sizeFactor_++] = kRow;
    2050                          } else {
    2051                               int kRow = iOriginalRow - numberTotal;
    2052                               CoinBigIndex start = rowStart[kRow];
    2053                               CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
    2054                               for (CoinBigIndex j = start; j < end; j++) {
    2055                                    int jRow = column[j];
    2056                                    int jNewRow = permuteInverse_[jRow];
    2057                                    if (jNewRow < iRow)
    2058                                         Arow[sizeFactor_++] = jNewRow;
    2059                               }
    2060                               // slack - should it be permute
    2061                               kRow = permuteInverse_[kRow+numberColumns];
    2062                               if (kRow < iRow)
    2063                                    Arow[sizeFactor_++] = kRow;
    2064                          }
    2065                          // Sort
    2066                          std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
    2067                     }
    2068                } else {
    2069                     // quadratic
    2070                     // transpose
    2071                     CoinPackedMatrix quadraticT;
    2072                     quadraticT.reverseOrderedCopyOf(*quadratic);
    2073                     const int * columnQuadratic = quadraticT.getIndices();
    2074                     const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
    2075                     const int * columnQuadraticLength = quadraticT.getVectorLengths();
    2076                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2077                          Astart[iRow] = sizeFactor_;
    2078                          int iOriginalRow = permute_[iRow];
    2079                          if (iOriginalRow < numberColumns) {
    2080                               // Quadratic bit
    2081                               CoinBigIndex j;
    2082                               for ( j = columnQuadraticStart[iOriginalRow];
    2083                                         j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
    2084                                    int jColumn = columnQuadratic[j];
    2085                                    int jNewColumn = permuteInverse_[jColumn];
    2086                                    if (jNewColumn < iRow)
    2087                                         Arow[sizeFactor_++] = jNewColumn;
    2088                               }
    2089                               // A may be upper triangular by mistake
    2090                               iColumn = iOriginalRow;
    2091                               CoinBigIndex start = columnStart[iColumn];
    2092                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2093                               for (j = start; j < end; j++) {
    2094                                    int kRow = row[j] + numberTotal;
    2095                                    kRow = permuteInverse_[kRow];
    2096                                    if (kRow < iRow)
    2097                                         Arow[sizeFactor_++] = kRow;
    2098                               }
    2099                          } else if (iOriginalRow < numberTotal) {
    2100                               int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
    2101                               if (kRow < iRow)
    2102                                    Arow[sizeFactor_++] = kRow;
    2103                          } else {
    2104                               int kRow = iOriginalRow - numberTotal;
    2105                               CoinBigIndex start = rowStart[kRow];
    2106                               CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
    2107                               for (CoinBigIndex j = start; j < end; j++) {
    2108                                    int jRow = column[j];
    2109                                    int jNewRow = permuteInverse_[jRow];
    2110                                    if (jNewRow < iRow)
    2111                                         Arow[sizeFactor_++] = jNewRow;
    2112                               }
    2113                               // slack - should it be permute
    2114                               kRow = permuteInverse_[kRow+numberColumns];
    2115                               if (kRow < iRow)
    2116                                    Arow[sizeFactor_++] = kRow;
    2117                          }
    2118                          // Sort
    2119                          std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
    2120                     }
    2121                }
    2122           } else {
    2123                if (!quadratic) {
    2124                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2125                          Astart[iRow] = sizeFactor_;
    2126                     }
    2127                } else {
    2128                     // Quadratic
    2129                     // transpose
    2130                     CoinPackedMatrix quadraticT;
    2131                     quadraticT.reverseOrderedCopyOf(*quadratic);
    2132                     const int * columnQuadratic = quadraticT.getIndices();
    2133                     const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
    2134                     const int * columnQuadraticLength = quadraticT.getVectorLengths();
    2135                     //const double * quadraticElement = quadraticT.getElements();
    2136                     for (iRow = 0; iRow < numberColumns; iRow++) {
    2137                          int iOriginalRow = permute_[iRow];
    2138                          Astart[iRow] = sizeFactor_;
    2139                          for (CoinBigIndex j = columnQuadraticStart[iOriginalRow];
    2140                                    j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
    2141                               int jColumn = columnQuadratic[j];
    2142                               int jNewColumn = permuteInverse_[jColumn];
    2143                               if (jNewColumn < iRow)
    2144                                    Arow[sizeFactor_++] = jNewColumn;
    2145                          }
    2146                     }
    2147                }
    2148                int iRow;
    2149                // slacks
    2150                for (iRow = 0; iRow < numberRowsModel; iRow++) {
    2151                     Astart[iRow+numberColumns] = sizeFactor_;
    2152                }
    2153                for (iRow = 0; iRow < numberRowsModel; iRow++) {
    2154                     Astart[iRow+numberTotal] = sizeFactor_;
    2155                     CoinBigIndex start = rowStart[iRow];
    2156                     CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
    2157                     for (CoinBigIndex j = start; j < end; j++) {
    2158                          Arow[sizeFactor_++] = column[j];
    2159                     }
    2160                     // slack
    2161                     Arow[sizeFactor_++] = numberColumns + iRow;
    2162                }
    2163           }
    2164           delete rowCopy;
    2165      }
    2166      Astart[numberRows_] = sizeFactor_;
    2167      firstDense_ = numberRows_;
    2168      symbolic1(Astart, Arow);
    2169      // Now fill in indices
    2170      try {
    2171           // too big
    2172           choleskyRow_ = new int[sizeFactor_];
    2173      } catch (...) {
    2174           // no memory
    2175           noMemory = true;
    2176      }
    2177      double sizeFactor = sizeFactor_;
    2178      if (!noMemory) {
    2179           // Do lower triangular
    2180           sizeFactor_ = 0;
    2181           int * which = Arow;
    2182           if (!doKKT_) {
    2183                for (iRow = 0; iRow < numberRows_; iRow++) {
    2184                     int number = 0;
    2185                     int iOriginalRow = permute_[iRow];
    2186                     Astart[iRow] = sizeFactor_;
    2187                     if (!rowsDropped_[iOriginalRow]) {
    2188                          CoinBigIndex startRow = rowStart[iOriginalRow];
    2189                          CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
    2190                          for (CoinBigIndex k = startRow; k < endRow; k++) {
    2191                               int iColumn = column[k];
    2192                               if (!whichDense_ || !whichDense_[iColumn]) {
    2193                                    CoinBigIndex start = columnStart[iColumn];
    2194                                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2195                                    for (CoinBigIndex j = start; j < end; j++) {
    2196                                         int jRow = row[j];
    2197                                         int jNewRow = permuteInverse_[jRow];
    2198                                         if (jNewRow > iRow && !rowsDropped_[jRow]) {
    2199                                              if (!used[jNewRow]) {
    2200                                                   used[jNewRow] = 1;
    2201                                                   which[number++] = jNewRow;
    2202                                              }
    2203                                         }
    2204                                    }
    2205                               }
    2206                          }
    2207                          sizeFactor_ += number;
    2208                          int j;
    2209                          for (j = 0; j < number; j++)
    2210                               used[which[j]] = 0;
    2211                          // Sort
    2212                          std::sort(which, which + number);
    2213                          // move which on
    2214                          which += number;
    2215                     }
    2216                }
    2217           } else {
    2218                // KKT
    2219                // temp
    2220                bool permuted = false;
    2221                for (iRow = 0; iRow < numberRows_; iRow++) {
    2222                     if (permute_[iRow] != iRow) {
    2223                          permuted = true;
    2224                          break;
    2225                     }
    2226                }
    2227                // but fake it
    2228                for (iRow = 0; iRow < numberRows_; iRow++) {
    2229                     //permute_[iRow]=iRow; // force no permute
    2230                     //permuteInverse_[permute_[iRow]]=iRow;
    2231                }
    2232                if (permuted) {
    2233                     // Need to permute - ugly
    2234                     if (!quadratic) {
    2235                          for (iRow = 0; iRow < numberRows_; iRow++) {
    2236                               Astart[iRow] = sizeFactor_;
    2237                               int iOriginalRow = permute_[iRow];
    2238                               if (iOriginalRow < numberColumns) {
    2239                                    iColumn = iOriginalRow;
    2240                                    CoinBigIndex start = columnStart[iColumn];
    2241                                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2242                                    for (CoinBigIndex j = start; j < end; j++) {
    2243                                         int kRow = row[j] + numberTotal;
    2244                                         kRow = permuteInverse_[kRow];
    2245                                         if (kRow > iRow)
    2246                                              Arow[sizeFactor_++] = kRow;
    2247                                    }
    2248                               } else if (iOriginalRow < numberTotal) {
    2249                                    int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
    2250                                    if (kRow > iRow)
    2251                                         Arow[sizeFactor_++] = kRow;
    2252                               } else {
    2253                                    int kRow = iOriginalRow - numberTotal;
    2254                                    CoinBigIndex start = rowStart[kRow];
    2255                                    CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
    2256                                    for (CoinBigIndex j = start; j < end; j++) {
    2257                                         int jRow = column[j];
    2258                                         int jNewRow = permuteInverse_[jRow];
    2259                                         if (jNewRow > iRow)
    2260                                              Arow[sizeFactor_++] = jNewRow;
    2261                                    }
    2262                                    // slack - should it be permute
    2263                                    kRow = permuteInverse_[kRow+numberColumns];
    2264                                    if (kRow > iRow)
    2265                                         Arow[sizeFactor_++] = kRow;
    2266                               }
    2267                               // Sort
    2268                               std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
    2269                          }
    2270                     } else {
    2271                          // quadratic
    2272                          const int * columnQuadratic = quadratic->getIndices();
    2273                          const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    2274                          const int * columnQuadraticLength = quadratic->getVectorLengths();
    2275                          for (iRow = 0; iRow < numberRows_; iRow++) {
    2276                               Astart[iRow] = sizeFactor_;
    2277                               int iOriginalRow = permute_[iRow];
    2278                               if (iOriginalRow < numberColumns) {
    2279                                    // Quadratic bit
    2280                                    CoinBigIndex j;
    2281                                    for ( j = columnQuadraticStart[iOriginalRow];
    2282                                              j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
    2283                                         int jColumn = columnQuadratic[j];
    2284                                         int jNewColumn = permuteInverse_[jColumn];
    2285                                         if (jNewColumn > iRow)
    2286                                              Arow[sizeFactor_++] = jNewColumn;
    2287                                    }
    2288                                    iColumn = iOriginalRow;
    2289                                    CoinBigIndex start = columnStart[iColumn];
    2290                                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2291                                    for (j = start; j < end; j++) {
    2292                                         int kRow = row[j] + numberTotal;
    2293                                         kRow = permuteInverse_[kRow];
    2294                                         if (kRow > iRow)
    2295                                              Arow[sizeFactor_++] = kRow;
    2296                                    }
    2297                               } else if (iOriginalRow < numberTotal) {
    2298                                    int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
    2299                                    if (kRow > iRow)
    2300                                         Arow[sizeFactor_++] = kRow;
    2301                               } else {
    2302                                    int kRow = iOriginalRow - numberTotal;
    2303                                    CoinBigIndex start = rowStart[kRow];
    2304                                    CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
    2305                                    for (CoinBigIndex j = start; j < end; j++) {
    2306                                         int jRow = column[j];
    2307                                         int jNewRow = permuteInverse_[jRow];
    2308                                         if (jNewRow > iRow)
    2309                                              Arow[sizeFactor_++] = jNewRow;
    2310                                    }
    2311                                    // slack - should it be permute
    2312                                    kRow = permuteInverse_[kRow+numberColumns];
    2313                                    if (kRow > iRow)
    2314                                         Arow[sizeFactor_++] = kRow;
    2315                               }
    2316                               // Sort
    2317                               std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
    2318                          }
    2319                     }
    2320                } else {
    2321                     if (!quadratic) {
    2322                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2323                               Astart[iColumn] = sizeFactor_;
    2324                               CoinBigIndex start = columnStart[iColumn];
    2325                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2326                               for (CoinBigIndex j = start; j < end; j++) {
    2327                                    Arow[sizeFactor_++] = row[j] + numberTotal;
    2328                               }
    2329                          }
    2330                     } else {
    2331                          // Quadratic
    2332                          const int * columnQuadratic = quadratic->getIndices();
    2333                          const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    2334                          const int * columnQuadraticLength = quadratic->getVectorLengths();
    2335                          //const double * quadraticElement = quadratic->getElements();
    2336                          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2337                               Astart[iColumn] = sizeFactor_;
    2338                               CoinBigIndex j;
    2339                               for ( j = columnQuadraticStart[iColumn];
    2340                                         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    2341                                    int jColumn = columnQuadratic[j];
    2342                                    if (jColumn > iColumn)
    2343                                         Arow[sizeFactor_++] = jColumn;
    2344                               }
    2345                               CoinBigIndex start = columnStart[iColumn];
    2346                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2347                               for ( j = start; j < end; j++) {
    2348                                    Arow[sizeFactor_++] = row[j] + numberTotal;
    2349                               }
    2350                          }
    2351                     }
    2352                     // slacks
    2353                     for (iRow = 0; iRow < numberRowsModel; iRow++) {
    2354                          Astart[iRow+numberColumns] = sizeFactor_;
    2355                          Arow[sizeFactor_++] = iRow + numberTotal;
    2356                     }
    2357                     // Transpose - nonzero diagonal (may regularize)
    2358                     for (iRow = 0; iRow < numberRowsModel; iRow++) {
    2359                          Astart[iRow+numberTotal] = sizeFactor_;
    2360                     }
    2361                }
    2362                Astart[numberRows_] = sizeFactor_;
    2363           }
    2364           symbolic2(Astart, Arow);
    2365           if (sizeIndex_ < sizeFactor_) {
    2366                int * indices = NULL;
    2367                try {
    2368                     indices = new int[sizeIndex_];
    2369                } catch (...) {
    2370                     // no memory
    2371                     noMemory = true;
    2372                }
    2373                if (!noMemory)  {
    2374                     CoinMemcpyN(choleskyRow_, sizeIndex_, indices);
    2375                     delete [] choleskyRow_;
    2376                     choleskyRow_ = indices;
    2377                }
    2378           }
    2379      }
    2380      delete [] used;
    2381      // Use cholesky regions
    2382      delete [] Astart;
    2383      delete [] Arow;
    2384      double flops = 0.0;
    2385      for (iRow = 0; iRow < numberRows_; iRow++) {
    2386           int length = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    2387           flops += static_cast<double> (length) * (length + 2.0);
    2388      }
    2389      if (model_->messageHandler()->logLevel() > 0)
    2390           std::cout << sizeFactor << " elements in sparse Cholesky, flop count " << flops << std::endl;
    2391      try {
    2392           sparseFactor_ = new longDouble [sizeFactor_];
    2393 #if CLP_LONG_CHOLESKY!=1
    2394           workDouble_ = new longDouble[numberRows_];
     2313            Astart[iRow] = sizeFactor_;
     2314            int iOriginalRow = permute_[iRow];
     2315            if (iOriginalRow < numberColumns) {
     2316              // Quadratic bit
     2317              CoinBigIndex j;
     2318              for (j = columnQuadraticStart[iOriginalRow];
     2319                   j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
     2320                int jColumn = columnQuadratic[j];
     2321                int jNewColumn = permuteInverse_[jColumn];
     2322                if (jNewColumn > iRow)
     2323                  Arow[sizeFactor_++] = jNewColumn;
     2324              }
     2325              iColumn = iOriginalRow;
     2326              CoinBigIndex start = columnStart[iColumn];
     2327              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2328              for (j = start; j < end; j++) {
     2329                int kRow = row[j] + numberTotal;
     2330                kRow = permuteInverse_[kRow];
     2331                if (kRow > iRow)
     2332                  Arow[sizeFactor_++] = kRow;
     2333              }
     2334            } else if (iOriginalRow < numberTotal) {
     2335              int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
     2336              if (kRow > iRow)
     2337                Arow[sizeFactor_++] = kRow;
     2338            } else {
     2339              int kRow = iOriginalRow - numberTotal;
     2340              CoinBigIndex start = rowStart[kRow];
     2341              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
     2342              for (CoinBigIndex j = start; j < end; j++) {
     2343                int jRow = column[j];
     2344                int jNewRow = permuteInverse_[jRow];
     2345                if (jNewRow > iRow)
     2346                  Arow[sizeFactor_++] = jNewRow;
     2347              }
     2348              // slack - should it be permute
     2349              kRow = permuteInverse_[kRow + numberColumns];
     2350              if (kRow > iRow)
     2351                Arow[sizeFactor_++] = kRow;
     2352            }
     2353            // Sort
     2354            std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
     2355          }
     2356        }
     2357      } else {
     2358        if (!quadratic) {
     2359          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2360            Astart[iColumn] = sizeFactor_;
     2361            CoinBigIndex start = columnStart[iColumn];
     2362            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2363            for (CoinBigIndex j = start; j < end; j++) {
     2364              Arow[sizeFactor_++] = row[j] + numberTotal;
     2365            }
     2366          }
     2367        } else {
     2368          // Quadratic
     2369          const int *columnQuadratic = quadratic->getIndices();
     2370          const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
     2371          const int *columnQuadraticLength = quadratic->getVectorLengths();
     2372          //const double * quadraticElement = quadratic->getElements();
     2373          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2374            Astart[iColumn] = sizeFactor_;
     2375            CoinBigIndex j;
     2376            for (j = columnQuadraticStart[iColumn];
     2377                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
     2378              int jColumn = columnQuadratic[j];
     2379              if (jColumn > iColumn)
     2380                Arow[sizeFactor_++] = jColumn;
     2381            }
     2382            CoinBigIndex start = columnStart[iColumn];
     2383            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2384            for (j = start; j < end; j++) {
     2385              Arow[sizeFactor_++] = row[j] + numberTotal;
     2386            }
     2387          }
     2388        }
     2389        // slacks
     2390        for (iRow = 0; iRow < numberRowsModel; iRow++) {
     2391          Astart[iRow + numberColumns] = sizeFactor_;
     2392          Arow[sizeFactor_++] = iRow + numberTotal;
     2393        }
     2394        // Transpose - nonzero diagonal (may regularize)
     2395        for (iRow = 0; iRow < numberRowsModel; iRow++) {
     2396          Astart[iRow + numberTotal] = sizeFactor_;
     2397        }
     2398      }
     2399      Astart[numberRows_] = sizeFactor_;
     2400    }
     2401    symbolic2(Astart, Arow);
     2402    if (sizeIndex_ < sizeFactor_) {
     2403      int *indices = NULL;
     2404      try {
     2405        indices = new int[sizeIndex_];
     2406      } catch (...) {
     2407        // no memory
     2408        noMemory = true;
     2409      }
     2410      if (!noMemory) {
     2411        CoinMemcpyN(choleskyRow_, sizeIndex_, indices);
     2412        delete[] choleskyRow_;
     2413        choleskyRow_ = indices;
     2414      }
     2415    }
     2416  }
     2417  delete[] used;
     2418  // Use cholesky regions
     2419  delete[] Astart;
     2420  delete[] Arow;
     2421  double flops = 0.0;
     2422  for (iRow = 0; iRow < numberRows_; iRow++) {
     2423    int length = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
     2424    flops += static_cast< double >(length) * (length + 2.0);
     2425  }
     2426  if (model_->messageHandler()->logLevel() > 0)
     2427    std::cout << sizeFactor << " elements in sparse Cholesky, flop count " << flops << std::endl;
     2428  try {
     2429    sparseFactor_ = new longDouble[sizeFactor_];
     2430#if CLP_LONG_CHOLESKY != 1
     2431    workDouble_ = new longDouble[numberRows_];
    23952432#else
    2396           // actually long double
    2397           workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]);
    2398 #endif
    2399           diagonal_ = new longDouble[numberRows_];
    2400      } catch (...) {
    2401           // no memory
    2402           noMemory = true;
    2403      }
    2404      if (noMemory) {
    2405           delete [] choleskyRow_;
    2406           choleskyRow_ = NULL;
    2407           delete [] choleskyStart_;
    2408           choleskyStart_ = NULL;
    2409           delete [] permuteInverse_;
    2410           permuteInverse_ = NULL;
    2411           delete [] permute_;
    2412           permute_ = NULL;
    2413           delete [] choleskyStart_;
    2414           choleskyStart_ = NULL;
    2415           delete [] indexStart_;
    2416           indexStart_ = NULL;
    2417           delete [] link_;
    2418           link_ = NULL;
    2419           delete [] workInteger_;
    2420           workInteger_ = NULL;
    2421           delete [] sparseFactor_;
    2422           sparseFactor_ = NULL;
    2423           delete [] workDouble_;
    2424           workDouble_ = NULL;
    2425           delete [] diagonal_;
    2426           diagonal_ = NULL;
    2427           delete [] clique_;
    2428           clique_ = NULL;
    2429           return -1;
    2430      }
    2431      return 0;
     2433    // actually long double
     2434    workDouble_ = reinterpret_cast< double * >(new CoinWorkDouble[numberRows_]);
     2435#endif
     2436    diagonal_ = new longDouble[numberRows_];
     2437  } catch (...) {
     2438    // no memory
     2439    noMemory = true;
     2440  }
     2441  if (noMemory) {
     2442    delete[] choleskyRow_;
     2443    choleskyRow_ = NULL;
     2444    delete[] choleskyStart_;
     2445    choleskyStart_ = NULL;
     2446    delete[] permuteInverse_;
     2447    permuteInverse_ = NULL;
     2448    delete[] permute_;
     2449    permute_ = NULL;
     2450    delete[] choleskyStart_;
     2451    choleskyStart_ = NULL;
     2452    delete[] indexStart_;
     2453    indexStart_ = NULL;
     2454    delete[] link_;
     2455    link_ = NULL;
     2456    delete[] workInteger_;
     2457    workInteger_ = NULL;
     2458    delete[] sparseFactor_;
     2459    sparseFactor_ = NULL;
     2460    delete[] workDouble_;
     2461    workDouble_ = NULL;
     2462    delete[] diagonal_;
     2463    diagonal_ = NULL;
     2464    delete[] clique_;
     2465    clique_ = NULL;
     2466    return -1;
     2467  }
     2468  return 0;
    24322469}
    2433 int
    2434 ClpCholeskyBase::symbolic1(const int * Astart, const int * Arow)
     2470int ClpCholeskyBase::symbolic1(const int *Astart, const int *Arow)
    24352471{
    2436      int * marked = reinterpret_cast<int *> (workInteger_);
    2437      int iRow;
    2438      // may not need to do this here but makes debugging easier
    2439      for (iRow = 0; iRow < numberRows_; iRow++) {
    2440           marked[iRow] = -1;
    2441           link_[iRow] = -1;
    2442           choleskyStart_[iRow] = 0; // counts
    2443      }
    2444      for (iRow = 0; iRow < numberRows_; iRow++) {
    2445           marked[iRow] = iRow;
    2446           for (CoinBigIndex j = Astart[iRow]; j < Astart[iRow+1]; j++) {
    2447                int kRow = Arow[j];
    2448                while (marked[kRow] != iRow ) {
    2449                     if (link_[kRow] < 0 )
    2450                          link_[kRow] = iRow;
    2451                     choleskyStart_[kRow]++;
    2452                     marked[kRow] = iRow;
    2453                     kRow = link_[kRow];
    2454                }
    2455           }
    2456      }
    2457      sizeFactor_ = 0;
    2458      for (iRow = 0; iRow < numberRows_; iRow++) {
    2459           int number = choleskyStart_[iRow];
    2460           choleskyStart_[iRow] = sizeFactor_;
    2461           sizeFactor_ += number;
    2462      }
    2463      choleskyStart_[numberRows_] = sizeFactor_;
    2464      return sizeFactor_;;
     2472  int *marked = reinterpret_cast< int * >(workInteger_);
     2473  int iRow;
     2474  // may not need to do this here but makes debugging easier
     2475  for (iRow = 0; iRow < numberRows_; iRow++) {
     2476    marked[iRow] = -1;
     2477    link_[iRow] = -1;
     2478    choleskyStart_[iRow] = 0; // counts
     2479  }
     2480  for (iRow = 0; iRow < numberRows_; iRow++) {
     2481    marked[iRow] = iRow;
     2482    for (CoinBigIndex j = Astart[iRow]; j < Astart[iRow + 1]; j++) {
     2483      int kRow = Arow[j];
     2484      while (marked[kRow] != iRow) {
     2485        if (link_[kRow] < 0)
     2486          link_[kRow] = iRow;
     2487        choleskyStart_[kRow]++;
     2488        marked[kRow] = iRow;
     2489        kRow = link_[kRow];
     2490      }
     2491    }
     2492  }
     2493  sizeFactor_ = 0;
     2494  for (iRow = 0; iRow < numberRows_; iRow++) {
     2495    int number = choleskyStart_[iRow];
     2496    choleskyStart_[iRow] = sizeFactor_;
     2497    sizeFactor_ += number;
     2498  }
     2499  choleskyStart_[numberRows_] = sizeFactor_;
     2500  return sizeFactor_;
     2501  ;
    24652502}
    2466 void
    2467 ClpCholeskyBase::symbolic2(const int * Astart, const int * Arow)
     2503void ClpCholeskyBase::symbolic2(const int *Astart, const int *Arow)
    24682504{
    2469      int * mergeLink = clique_;
    2470      int * marker = reinterpret_cast<int *> (workInteger_);
    2471      int iRow;
    2472      for (iRow = 0; iRow < numberRows_; iRow++) {
    2473           marker[iRow] = -1;
    2474           mergeLink[iRow] = -1;
    2475           link_[iRow] = -1; // not needed but makes debugging easier
    2476      }
    2477      int start = 0;
    2478      int end = 0;
    2479      choleskyStart_[0] = 0;
     2505  int *mergeLink = clique_;
     2506  int *marker = reinterpret_cast< int * >(workInteger_);
     2507  int iRow;
     2508  for (iRow = 0; iRow < numberRows_; iRow++) {
     2509    marker[iRow] = -1;
     2510    mergeLink[iRow] = -1;
     2511    link_[iRow] = -1; // not needed but makes debugging easier
     2512  }
     2513  int start = 0;
     2514  int end = 0;
     2515  choleskyStart_[0] = 0;
    24802516
    2481      for (iRow = 0; iRow < numberRows_; iRow++) {
    2482           int nz = 0;
    2483           int merge = mergeLink[iRow];
    2484           bool marked = false;
    2485           if (merge < 0)
    2486                marker[iRow] = iRow;
    2487           else
    2488                marker[iRow] = merge;
    2489           start = end;
    2490           int startSub = start;
    2491           link_[iRow] = numberRows_;
    2492           CoinBigIndex j;
    2493           for ( j = Astart[iRow]; j < Astart[iRow+1]; j++) {
    2494                int kRow = Arow[j];
    2495                int k = iRow;
    2496                int linked = link_[iRow];
     2517  for (iRow = 0; iRow < numberRows_; iRow++) {
     2518    int nz = 0;
     2519    int merge = mergeLink[iRow];
     2520    bool marked = false;
     2521    if (merge < 0)
     2522      marker[iRow] = iRow;
     2523    else
     2524      marker[iRow] = merge;
     2525    start = end;
     2526    int startSub = start;
     2527    link_[iRow] = numberRows_;
     2528    CoinBigIndex j;
     2529    for (j = Astart[iRow]; j < Astart[iRow + 1]; j++) {
     2530      int kRow = Arow[j];
     2531      int k = iRow;
     2532      int linked = link_[iRow];
    24972533#ifndef NDEBUG
    2498                int count = 0;
    2499 #endif
    2500                while (linked <= kRow) {
    2501                     k = linked;
    2502                     linked = link_[k];
     2534      int count = 0;
     2535#endif
     2536      while (linked <= kRow) {
     2537        k = linked;
     2538        linked = link_[k];
    25032539#ifndef NDEBUG
    2504                     count++;
    2505                     assert (count < numberRows_);
    2506 #endif
    2507                }
    2508                nz++;
    2509                link_[k] = kRow;
    2510                link_[kRow] = linked;
    2511                if (marker[kRow] != marker[iRow])
    2512                     marked = true;
    2513           }
    2514           bool reuse = false;
    2515           // Check if we can re-use indices
    2516           if (!marked && merge >= 0 && mergeLink[merge] < 0) {
    2517                // can re-use all
    2518                startSub = indexStart_[merge] + 1;
    2519                nz = choleskyStart_[merge+1] - (choleskyStart_[merge] + 1);
    2520                reuse = true;
    2521           } else {
    2522                // See if we can re-use any
    2523                int k = mergeLink[iRow];
    2524                int maxLength = 0;
    2525                while (k >= 0) {
    2526                     int length = choleskyStart_[k+1] - (choleskyStart_[k] + 1);
    2527                     int start = indexStart_[k] + 1;
    2528                     int stop = start + length;
    2529                     if (length > maxLength) {
    2530                          maxLength = length;
    2531                          startSub = start;
    2532                     }
    2533                     int linked = iRow;
    2534                     for (CoinBigIndex j = start; j < stop; j++) {
    2535                          int kRow = choleskyRow_[j];
    2536                          int kk = linked;
    2537                          linked = link_[kk];
    2538                          while (linked < kRow) {
    2539                               kk = linked;
    2540                               linked = link_[kk];
    2541                          }
    2542                          if (linked != kRow) {
    2543                               nz++;
    2544                               link_[kk] = kRow;
    2545                               link_[kRow] = linked;
    2546                               linked = kRow;
    2547                          }
    2548                     }
    2549                     k = mergeLink[k];
    2550                }
    2551                if (nz == maxLength)
    2552                     reuse = true; // can re-use
    2553           }
    2554           //reuse=false; //temp
    2555           if (!reuse) {
    2556                end += nz;
    2557                startSub = start;
    2558                int kRow = iRow;
    2559                for (int j = start; j < end; j++) {
    2560                     kRow = link_[kRow];
    2561                     choleskyRow_[j] = kRow;
    2562                     assert (kRow < numberRows_);
    2563                     marker[kRow] = iRow;
    2564                }
    2565                marker[iRow] = iRow;
    2566           }
    2567           indexStart_[iRow] = startSub;
    2568           choleskyStart_[iRow+1] = choleskyStart_[iRow] + nz;
    2569           if (nz > 1) {
    2570                int kRow = choleskyRow_[startSub];
    2571                mergeLink[iRow] = mergeLink[kRow];
    2572                mergeLink[kRow] = iRow;
    2573           }
    2574           // should not be needed
    2575           //std::sort(choleskyRow_+indexStart_[iRow]
    2576           //      ,choleskyRow_+indexStart_[iRow]+nz);
    2577           //#define CLP_DEBUG
     2540        count++;
     2541        assert(count < numberRows_);
     2542#endif
     2543      }
     2544      nz++;
     2545      link_[k] = kRow;
     2546      link_[kRow] = linked;
     2547      if (marker[kRow] != marker[iRow])
     2548        marked = true;
     2549    }
     2550    bool reuse = false;
     2551    // Check if we can re-use indices
     2552    if (!marked && merge >= 0 && mergeLink[merge] < 0) {
     2553      // can re-use all
     2554      startSub = indexStart_[merge] + 1;
     2555      nz = choleskyStart_[merge + 1] - (choleskyStart_[merge] + 1);
     2556      reuse = true;
     2557    } else {
     2558      // See if we can re-use any
     2559      int k = mergeLink[iRow];
     2560      int maxLength = 0;
     2561      while (k >= 0) {
     2562        int length = choleskyStart_[k + 1] - (choleskyStart_[k] + 1);
     2563        int start = indexStart_[k] + 1;
     2564        int stop = start + length;
     2565        if (length > maxLength) {
     2566          maxLength = length;
     2567          startSub = start;
     2568        }
     2569        int linked = iRow;
     2570        for (CoinBigIndex j = start; j < stop; j++) {
     2571          int kRow = choleskyRow_[j];
     2572          int kk = linked;
     2573          linked = link_[kk];
     2574          while (linked < kRow) {
     2575            kk = linked;
     2576            linked = link_[kk];
     2577          }
     2578          if (linked != kRow) {
     2579            nz++;
     2580            link_[kk] = kRow;
     2581            link_[kRow] = linked;
     2582            linked = kRow;
     2583          }
     2584        }
     2585        k = mergeLink[k];
     2586      }
     2587      if (nz == maxLength)
     2588        reuse = true; // can re-use
     2589    }
     2590    //reuse=false; //temp
     2591    if (!reuse) {
     2592      end += nz;
     2593      startSub = start;
     2594      int kRow = iRow;
     2595      for (int j = start; j < end; j++) {
     2596        kRow = link_[kRow];
     2597        choleskyRow_[j] = kRow;
     2598        assert(kRow < numberRows_);
     2599        marker[kRow] = iRow;
     2600      }
     2601      marker[iRow] = iRow;
     2602    }
     2603    indexStart_[iRow] = startSub;
     2604    choleskyStart_[iRow + 1] = choleskyStart_[iRow] + nz;
     2605    if (nz > 1) {
     2606      int kRow = choleskyRow_[startSub];
     2607      mergeLink[iRow] = mergeLink[kRow];
     2608      mergeLink[kRow] = iRow;
     2609    }
     2610    // should not be needed
     2611    //std::sort(choleskyRow_+indexStart_[iRow]
     2612    //      ,choleskyRow_+indexStart_[iRow]+nz);
     2613    //#define CLP_DEBUG
    25782614#ifdef CLP_DEBUG
    2579           int last = -1;
    2580           for ( j = indexStart_[iRow]; j < indexStart_[iRow] + nz; j++) {
    2581                int kRow = choleskyRow_[j];
    2582                assert (kRow > last);
    2583                last = kRow;
    2584           }
    2585 #endif
    2586      }
    2587      sizeFactor_ = choleskyStart_[numberRows_];
    2588      sizeIndex_ = start;
    2589      // find dense segment here
    2590      int numberleft = numberRows_;
    2591      for (iRow = 0; iRow < numberRows_; iRow++) {
    2592           CoinBigIndex left = sizeFactor_ - choleskyStart_[iRow];
    2593           double n = numberleft;
    2594           double threshold = n * (n - 1.0) * 0.5 * goDense_;
    2595           if ( left >= threshold)
    2596                break;
    2597           numberleft--;
    2598      }
    2599      //iRow=numberRows_;
    2600      int nDense = numberRows_ - iRow;
     2615    int last = -1;
     2616    for (j = indexStart_[iRow]; j < indexStart_[iRow] + nz; j++) {
     2617      int kRow = choleskyRow_[j];
     2618      assert(kRow > last);
     2619      last = kRow;
     2620    }
     2621#endif
     2622  }
     2623  sizeFactor_ = choleskyStart_[numberRows_];
     2624  sizeIndex_ = start;
     2625  // find dense segment here
     2626  int numberleft = numberRows_;
     2627  for (iRow = 0; iRow < numberRows_; iRow++) {
     2628    CoinBigIndex left = sizeFactor_ - choleskyStart_[iRow];
     2629    double n = numberleft;
     2630    double threshold = n * (n - 1.0) * 0.5 * goDense_;
     2631    if (left >= threshold)
     2632      break;
     2633    numberleft--;
     2634  }
     2635  //iRow=numberRows_;
     2636  int nDense = numberRows_ - iRow;
    26012637#define DENSE_THRESHOLD 8
    2602      // don't do if dense columns
    2603      if (nDense >= DENSE_THRESHOLD && !dense_) {
    2604        COIN_DETAIL_PRINT(printf("Going dense for last %d rows\n", nDense));
    2605           // make sure we don't disturb any indices
    2606           int k = 0;
    2607           for (int jRow = 0; jRow < iRow; jRow++) {
    2608                int nz = choleskyStart_[jRow+1] - choleskyStart_[jRow];
    2609                k = CoinMax(k, indexStart_[jRow] + nz);
    2610           }
    2611           indexStart_[iRow] = k;
    2612           int j;
    2613           for (j = iRow + 1; j < numberRows_; j++) {
    2614                choleskyRow_[k++] = j;
    2615                indexStart_[j] = k;
    2616           }
    2617           sizeIndex_ = k;
    2618           assert (k <= sizeFactor_); // can't happen with any reasonable defaults
    2619           k = choleskyStart_[iRow];
    2620           for (j = iRow + 1; j <= numberRows_; j++) {
    2621                k += numberRows_ - j;
    2622                choleskyStart_[j] = k;
    2623           }
    2624           // allow for blocked dense
    2625           ClpCholeskyDense dense;
    2626           sizeFactor_ = choleskyStart_[iRow] + dense.space(nDense);
    2627           firstDense_ = iRow;
    2628           if (doKKT_) {
    2629                // redo permute so negative ones first
    2630                int putN = firstDense_;
    2631                int putP = 0;
    2632                int numberRowsModel = model_->numberRows();
    2633                int numberColumns = model_->numberColumns();
    2634                int numberTotal = numberColumns + numberRowsModel;
    2635                for (iRow = firstDense_; iRow < numberRows_; iRow++) {
    2636                     int originalRow = permute_[iRow];
    2637                     if (originalRow < numberTotal)
    2638                          permute_[putN++] = originalRow;
    2639                     else
    2640                          permuteInverse_[putP++] = originalRow;
    2641                }
    2642                for (iRow = putN; iRow < numberRows_; iRow++) {
    2643                     permute_[iRow] = permuteInverse_[iRow-putN];
    2644                }
    2645                for (iRow = 0; iRow < numberRows_; iRow++) {
    2646                     permuteInverse_[permute_[iRow]] = iRow;
    2647                }
    2648           }
    2649      }
    2650      // Clean up clique info
    2651      for (iRow = 0; iRow < numberRows_; iRow++)
    2652           clique_[iRow] = 0;
    2653      int lastClique = -1;
    2654      bool inClique = false;
    2655      for (iRow = 1; iRow < firstDense_; iRow++) {
    2656           int sizeLast = choleskyStart_[iRow] - choleskyStart_[iRow-1];
    2657           int sizeThis = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    2658           if (indexStart_[iRow] == indexStart_[iRow-1] + 1 &&
    2659                     sizeThis == sizeLast - 1 &&
    2660                     sizeThis) {
    2661                // in clique
    2662                if (!inClique) {
    2663                     inClique = true;
    2664                     lastClique = iRow - 1;
    2665                }
    2666           } else if (inClique) {
    2667                int sizeClique = iRow - lastClique;
    2668                for (int i = lastClique; i < iRow; i++) {
    2669                     clique_[i] = sizeClique;
    2670                     sizeClique--;
    2671                }
    2672                inClique = false;
    2673           }
    2674      }
    2675      if (inClique) {
    2676           int sizeClique = iRow - lastClique;
    2677           for (int i = lastClique; i < iRow; i++) {
    2678                clique_[i] = sizeClique;
    2679                sizeClique--;
    2680           }
    2681      }
    2682      //for (iRow=0;iRow<numberRows_;iRow++)
    2683      //clique_[iRow]=0;
     2638  // don't do if dense columns
     2639  if (nDense >= DENSE_THRESHOLD && !dense_) {
     2640    COIN_DETAIL_PRINT(printf("Going dense for last %d rows\n", nDense));
     2641    // make sure we don't disturb any indices
     2642    int k = 0;
     2643    for (int jRow = 0; jRow < iRow; jRow++) {
     2644      int nz = choleskyStart_[jRow + 1] - choleskyStart_[jRow];
     2645      k = CoinMax(k, indexStart_[jRow] + nz);
     2646    }
     2647    indexStart_[iRow] = k;
     2648    int j;
     2649    for (j = iRow + 1; j < numberRows_; j++) {
     2650      choleskyRow_[k++] = j;
     2651      indexStart_[j] = k;
     2652    }
     2653    sizeIndex_ = k;
     2654    assert(k <= sizeFactor_); // can't happen with any reasonable defaults
     2655    k = choleskyStart_[iRow];
     2656    for (j = iRow + 1; j <= numberRows_; j++) {
     2657      k += numberRows_ - j;
     2658      choleskyStart_[j] = k;
     2659    }
     2660    // allow for blocked dense
     2661    ClpCholeskyDense dense;
     2662    sizeFactor_ = choleskyStart_[iRow] + dense.space(nDense);
     2663    firstDense_ = iRow;
     2664    if (doKKT_) {
     2665      // redo permute so negative ones first
     2666      int putN = firstDense_;
     2667      int putP = 0;
     2668      int numberRowsModel = model_->numberRows();
     2669      int numberColumns = model_->numberColumns();
     2670      int numberTotal = numberColumns + numberRowsModel;
     2671      for (iRow = firstDense_; iRow < numberRows_; iRow++) {
     2672        int originalRow = permute_[iRow];
     2673        if (originalRow < numberTotal)
     2674          permute_[putN++] = originalRow;
     2675        else
     2676          permuteInverse_[putP++] = originalRow;
     2677      }
     2678      for (iRow = putN; iRow < numberRows_; iRow++) {
     2679        permute_[iRow] = permuteInverse_[iRow - putN];
     2680      }
     2681      for (iRow = 0; iRow < numberRows_; iRow++) {
     2682        permuteInverse_[permute_[iRow]] = iRow;
     2683      }
     2684    }
     2685  }
     2686  // Clean up clique info
     2687  for (iRow = 0; iRow < numberRows_; iRow++)
     2688    clique_[iRow] = 0;
     2689  int lastClique = -1;
     2690  bool inClique = false;
     2691  for (iRow = 1; iRow < firstDense_; iRow++) {
     2692    int sizeLast = choleskyStart_[iRow] - choleskyStart_[iRow - 1];
     2693    int sizeThis = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
     2694    if (indexStart_[iRow] == indexStart_[iRow - 1] + 1 && sizeThis == sizeLast - 1 && sizeThis) {
     2695      // in clique
     2696      if (!inClique) {
     2697        inClique = true;
     2698        lastClique = iRow - 1;
     2699      }
     2700    } else if (inClique) {
     2701      int sizeClique = iRow - lastClique;
     2702      for (int i = lastClique; i < iRow; i++) {
     2703        clique_[i] = sizeClique;
     2704        sizeClique--;
     2705      }
     2706      inClique = false;
     2707    }
     2708  }
     2709  if (inClique) {
     2710    int sizeClique = iRow - lastClique;
     2711    for (int i = lastClique; i < iRow; i++) {
     2712      clique_[i] = sizeClique;
     2713      sizeClique--;
     2714    }
     2715  }
     2716  //for (iRow=0;iRow<numberRows_;iRow++)
     2717  //clique_[iRow]=0;
    26842718}
    26852719/* Factorize - filling in rowsDropped and returning number dropped */
    2686 int
    2687 ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped)
     2720int ClpCholeskyBase::factorize(const CoinWorkDouble *diagonal, int *rowsDropped)
    26882721{
    2689      const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    2690      const int * columnLength = model_->clpMatrix()->getVectorLengths();
    2691      const int * row = model_->clpMatrix()->getIndices();
    2692      const double * element = model_->clpMatrix()->getElements();
    2693      const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    2694      const int * rowLength = rowCopy_->getVectorLengths();
    2695      const int * column = rowCopy_->getIndices();
    2696      const double * elementByRow = rowCopy_->getElements();
    2697      int numberColumns = model_->clpMatrix()->getNumCols();
    2698      //perturbation
    2699      CoinWorkDouble perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
    2700      //perturbation=perturbation*perturbation*100000000.0;
    2701      if (perturbation > 1.0) {
     2722  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
     2723  const int *columnLength = model_->clpMatrix()->getVectorLengths();
     2724  const int *row = model_->clpMatrix()->getIndices();
     2725  const double *element = model_->clpMatrix()->getElements();
     2726  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
     2727  const int *rowLength = rowCopy_->getVectorLengths();
     2728  const int *column = rowCopy_->getIndices();
     2729  const double *elementByRow = rowCopy_->getElements();
     2730  int numberColumns = model_->clpMatrix()->getNumCols();
     2731  //perturbation
     2732  CoinWorkDouble perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
     2733  //perturbation=perturbation*perturbation*100000000.0;
     2734  if (perturbation > 1.0) {
    27022735#ifdef COIN_DEVELOP
    2703           //if (model_->model()->logLevel()&4)
    2704           std::cout << "large perturbation " << perturbation << std::endl;
    2705 #endif
    2706           perturbation = CoinSqrt(perturbation);
    2707           perturbation = 1.0;
    2708      }
    2709      int iRow;
    2710      int iColumn;
    2711      longDouble * work = workDouble_;
    2712      CoinZeroN(work, numberRows_);
    2713      int newDropped = 0;
    2714      CoinWorkDouble largest = 1.0;
    2715      CoinWorkDouble smallest = COIN_DBL_MAX;
    2716      int numberDense = 0;
    2717      if (!doKKT_) {
    2718           const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
    2719           if (dense_)
    2720                numberDense = dense_->numberRows();
    2721           if (whichDense_) {
    2722                longDouble * denseDiagonal = dense_->diagonal();
    2723                longDouble * dense = denseColumn_;
    2724                int iDense = 0;
    2725                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    2726                     if (whichDense_[iColumn]) {
    2727                          CoinZeroN(dense, numberRows_);
    2728                          CoinBigIndex start = columnStart[iColumn];
    2729                          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2730                          if (diagonal[iColumn]) {
    2731                               denseDiagonal[iDense++] = 1.0 / diagonal[iColumn];
    2732                               for (CoinBigIndex j = start; j < end; j++) {
    2733                                    int jRow = row[j];
    2734                                    dense[jRow] = element[j];
    2735                               }
    2736                          } else {
    2737                               denseDiagonal[iDense++] = 1.0;
    2738                          }
    2739                          dense += numberRows_;
    2740                     }
    2741                }
    2742           }
    2743           CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
    2744           delta2 *= delta2;
    2745           // largest in initial matrix
    2746           CoinWorkDouble largest2 = 1.0e-20;
    2747           for (iRow = 0; iRow < numberRows_; iRow++) {
    2748                longDouble * put = sparseFactor_ + choleskyStart_[iRow];
    2749                int * which = choleskyRow_ + indexStart_[iRow];
    2750                int iOriginalRow = permute_[iRow];
    2751                int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    2752                if (!rowLength[iOriginalRow])
    2753                     rowsDropped_[iOriginalRow] = 1;
    2754                if (!rowsDropped_[iOriginalRow]) {
    2755                     CoinBigIndex startRow = rowStart[iOriginalRow];
    2756                     CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
    2757                     work[iRow] = diagonalSlack[iOriginalRow] + delta2;
    2758                     for (CoinBigIndex k = startRow; k < endRow; k++) {
    2759                          int iColumn = column[k];
    2760                          if (!whichDense_ || !whichDense_[iColumn]) {
    2761                               CoinBigIndex start = columnStart[iColumn];
    2762                               CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2763                               CoinWorkDouble multiplier = diagonal[iColumn] * elementByRow[k];
    2764                               for (CoinBigIndex j = start; j < end; j++) {
    2765                                    int jRow = row[j];
    2766                                    int jNewRow = permuteInverse_[jRow];
    2767                                    if (jNewRow >= iRow && !rowsDropped_[jRow]) {
    2768                                         CoinWorkDouble value = element[j] * multiplier;
    2769                                         work[jNewRow] += value;
    2770                                    }
    2771                               }
    2772                          }
    2773                     }
    2774                     diagonal_[iRow] = work[iRow];
    2775                     largest2 = CoinMax(largest2, CoinAbs(work[iRow]));
    2776                     work[iRow] = 0.0;
    2777                     int j;
    2778                     for (j = 0; j < number; j++) {
    2779                          int jRow = which[j];
    2780                          put[j] = work[jRow];
    2781                          largest2 = CoinMax(largest2, CoinAbs(work[jRow]));
    2782                          work[jRow] = 0.0;
    2783                     }
    2784                } else {
    2785                     // dropped
    2786                     diagonal_[iRow] = 1.0;
    2787                     int j;
    2788                     for (j = 1; j < number; j++) {
    2789                          put[j] = 0.0;
    2790                     }
    2791                }
    2792           }
    2793           //check sizes
    2794           largest2 *= 1.0e-20;
    2795           largest = CoinMin(largest2, CHOL_SMALL_VALUE);
    2796           int numberDroppedBefore = 0;
    2797           for (iRow = 0; iRow < numberRows_; iRow++) {
    2798                int dropped = rowsDropped_[iRow];
    2799                // Move to int array
    2800                rowsDropped[iRow] = dropped;
    2801                if (!dropped) {
    2802                     CoinWorkDouble diagonal = diagonal_[iRow];
    2803                     if (diagonal > largest2) {
    2804                          diagonal_[iRow] = diagonal + perturbation;
    2805                     } else {
    2806                          diagonal_[iRow] = diagonal + perturbation;
    2807                          rowsDropped[iRow] = 2;
    2808                          numberDroppedBefore++;
    2809                          //printf("dropped - small diagonal %g\n",diagonal);
    2810                     }
    2811                }
    2812           }
    2813           doubleParameters_[10] = CoinMax(1.0e-20, largest);
    2814           integerParameters_[20] = 0;
    2815           doubleParameters_[3] = 0.0;
    2816           doubleParameters_[4] = COIN_DBL_MAX;
    2817           integerParameters_[34] = 0; // say all must be positive
    2818           factorizePart2(rowsDropped);
    2819           newDropped = integerParameters_[20] + numberDroppedBefore;
    2820           largest = doubleParameters_[3];
    2821           smallest = doubleParameters_[4];
    2822           if (model_->messageHandler()->logLevel() > 1)
    2823                std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
    2824           choleskyCondition_ = largest / smallest;
    2825           if (whichDense_) {
    2826                int i;
    2827                for ( i = 0; i < numberRows_; i++) {
    2828                     assert (diagonal_[i] >= 0.0);
    2829                     diagonal_[i] = CoinSqrt(diagonal_[i]);
    2830                }
    2831                // Update dense columns (just L)
    2832                // Zero out dropped rows
    2833                for (i = 0; i < numberDense; i++) {
    2834                     longDouble * a = denseColumn_ + i * numberRows_;
    2835                     for (int j = 0; j < numberRows_; j++) {
    2836                          if (rowsDropped[j])
    2837                               a[j] = 0.0;
    2838                     }
    2839                     for (i = 0; i < numberRows_; i++) {
    2840                          int iRow = permute_[i];
    2841                          workDouble_[i] = a[iRow];
    2842                     }
    2843                     for (i = 0; i < numberRows_; i++) {
    2844                          CoinWorkDouble value = workDouble_[i];
    2845                          CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
    2846                          CoinBigIndex j;
    2847                          for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
    2848                               int iRow = choleskyRow_[j+offset];
    2849                               workDouble_[iRow] -= sparseFactor_[j] * value;
    2850                          }
    2851                     }
    2852                     for (i = 0; i < numberRows_; i++) {
    2853                          int iRow = permute_[i];
    2854                          a[iRow] = workDouble_[i] * diagonal_[i];
    2855                     }
    2856                }
    2857                dense_->resetRowsDropped();
    2858                longDouble * denseBlob = dense_->aMatrix();
    2859                longDouble * denseDiagonal = dense_->diagonal();
    2860                // Update dense matrix
    2861                for (i = 0; i < numberDense; i++) {
    2862                     const longDouble * a = denseColumn_ + i * numberRows_;
    2863                     // do diagonal
    2864                     CoinWorkDouble value = denseDiagonal[i];
    2865                     const longDouble * b = denseColumn_ + i * numberRows_;
    2866                     for (int k = 0; k < numberRows_; k++)
    2867                          value += a[k] * b[k];
    2868                     denseDiagonal[i] = value;
    2869                     for (int j = i + 1; j < numberDense; j++) {
    2870                          CoinWorkDouble value = 0.0;
    2871                          const longDouble * b = denseColumn_ + j * numberRows_;
    2872                          for (int k = 0; k < numberRows_; k++)
    2873                               value += a[k] * b[k];
    2874                          *denseBlob = value;
    2875                          denseBlob++;
    2876                     }
    2877                }
    2878                // dense cholesky (? long double)
    2879                int * dropped = new int [numberDense];
    2880                dense_->factorizePart2(dropped);
    2881                delete [] dropped;
    2882           }
    2883           // try allowing all every time
    2884           //printf("trying ?\n");
    2885           //for (iRow=0;iRow<numberRows_;iRow++) {
    2886           //rowsDropped[iRow]=0;
    2887           //rowsDropped_[iRow]=0;
    2888           //}
    2889           bool cleanCholesky;
    2890           //if (model_->numberIterations()<20||(model_->numberIterations()&1)==0)
    2891           if (model_->numberIterations() < 2000)
    2892                cleanCholesky = true;
    2893           else
    2894                cleanCholesky = false;
    2895           if (cleanCholesky) {
    2896                //drop fresh makes some formADAT easier
    2897                if (newDropped || numberRowsDropped_) {
    2898                     newDropped = 0;
    2899                     for (int i = 0; i < numberRows_; i++) {
    2900                          char dropped = static_cast<char>(rowsDropped[i]);
    2901                          rowsDropped_[i] = dropped;
    2902                          rowsDropped_[i] = 0;
    2903                          if (dropped == 2) {
    2904                               //dropped this time
    2905                               rowsDropped[newDropped++] = i;
    2906                               rowsDropped_[i] = 0;
    2907                          }
    2908                     }
    2909                     numberRowsDropped_ = newDropped;
    2910                     newDropped = -(2 + newDropped);
    2911                }
     2736    //if (model_->model()->logLevel()&4)
     2737    std::cout << "large perturbation " << perturbation << std::endl;
     2738#endif
     2739    perturbation = CoinSqrt(perturbation);
     2740    perturbation = 1.0;
     2741  }
     2742  int iRow;
     2743  int iColumn;
     2744  longDouble *work = workDouble_;
     2745  CoinZeroN(work, numberRows_);
     2746  int newDropped = 0;
     2747  CoinWorkDouble largest = 1.0;
     2748  CoinWorkDouble smallest = COIN_DBL_MAX;
     2749  int numberDense = 0;
     2750  if (!doKKT_) {
     2751    const CoinWorkDouble *diagonalSlack = diagonal + numberColumns;
     2752    if (dense_)
     2753      numberDense = dense_->numberRows();
     2754    if (whichDense_) {
     2755      longDouble *denseDiagonal = dense_->diagonal();
     2756      longDouble *dense = denseColumn_;
     2757      int iDense = 0;
     2758      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     2759        if (whichDense_[iColumn]) {
     2760          CoinZeroN(dense, numberRows_);
     2761          CoinBigIndex start = columnStart[iColumn];
     2762          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2763          if (diagonal[iColumn]) {
     2764            denseDiagonal[iDense++] = 1.0 / diagonal[iColumn];
     2765            for (CoinBigIndex j = start; j < end; j++) {
     2766              int jRow = row[j];
     2767              dense[jRow] = element[j];
     2768            }
    29122769          } else {
    2913                if (newDropped) {
    2914                     newDropped = 0;
    2915                     for (int i = 0; i < numberRows_; i++) {
    2916                          char dropped = static_cast<char>(rowsDropped[i]);
    2917                          rowsDropped_[i] = dropped;
    2918                          if (dropped == 2) {
    2919                               //dropped this time
    2920                               rowsDropped[newDropped++] = i;
    2921                               rowsDropped_[i] = 1;
    2922                          }
    2923                     }
    2924                }
    2925                numberRowsDropped_ += newDropped;
    2926                if (numberRowsDropped_ && 0) {
    2927                     std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
    2928                               numberRowsDropped_ << " dropped)";
    2929                     if (newDropped) {
    2930                          std::cout << " ( " << newDropped << " dropped this time)";
    2931                     }
    2932                     std::cout << std::endl;
    2933                }
    2934           }
    2935      } else {
    2936           //KKT
    2937           CoinPackedMatrix * quadratic = NULL;
    2938           ClpQuadraticObjective * quadraticObj =
    2939                (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
    2940           if (quadraticObj)
    2941                quadratic = quadraticObj->quadraticObjective();
    2942           // matrix
    2943           int numberRowsModel = model_->numberRows();
    2944           int numberColumns = model_->numberColumns();
    2945           int numberTotal = numberColumns + numberRowsModel;
    2946           // temp
    2947           bool permuted = false;
    2948           for (iRow = 0; iRow < numberRows_; iRow++) {
    2949                if (permute_[iRow] != iRow) {
    2950                     permuted = true;
    2951                     break;
    2952                }
    2953           }
    2954           // but fake it
    2955           for (iRow = 0; iRow < numberRows_; iRow++) {
    2956                //permute_[iRow]=iRow; // force no permute
    2957                //permuteInverse_[permute_[iRow]]=iRow;
    2958           }
    2959           if (permuted) {
    2960                CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
    2961                delta2 *= delta2;
    2962                // Need to permute - ugly
    2963                if (!quadratic) {
    2964                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2965                          longDouble * put = sparseFactor_ + choleskyStart_[iRow];
    2966                          int * which = choleskyRow_ + indexStart_[iRow];
    2967                          int iOriginalRow = permute_[iRow];
    2968                          if (iOriginalRow < numberColumns) {
    2969                               iColumn = iOriginalRow;
    2970                               CoinWorkDouble value = diagonal[iColumn];
    2971                               if (CoinAbs(value) > 1.0e-100) {
    2972                                    value = 1.0 / value;
    2973                                    largest = CoinMax(largest, CoinAbs(value));
    2974                                    diagonal_[iRow] = -value;
    2975                                    CoinBigIndex start = columnStart[iColumn];
    2976                                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    2977                                    for (CoinBigIndex j = start; j < end; j++) {
    2978                                         int kRow = row[j] + numberTotal;
    2979                                         kRow = permuteInverse_[kRow];
    2980                                         if (kRow > iRow) {
    2981                                              work[kRow] = element[j];
    2982                                              largest = CoinMax(largest, CoinAbs(element[j]));
    2983                                         }
    2984                                    }
    2985                               } else {
    2986                                    diagonal_[iRow] = -value;
    2987                               }
    2988                          } else if (iOriginalRow < numberTotal) {
    2989                               CoinWorkDouble value = diagonal[iOriginalRow];
    2990                               if (CoinAbs(value) > 1.0e-100) {
    2991                                    value = 1.0 / value;
    2992                                    largest = CoinMax(largest, CoinAbs(value));
    2993                               } else {
    2994                                    value = 1.0e100;
    2995                               }
    2996                               diagonal_[iRow] = -value;
    2997                               int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
    2998                               if (kRow > iRow)
    2999                                    work[kRow] = -1.0;
    3000                          } else {
    3001                               diagonal_[iRow] = delta2;
    3002                               int kRow = iOriginalRow - numberTotal;
    3003                               CoinBigIndex start = rowStart[kRow];
    3004                               CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
    3005                               for (CoinBigIndex j = start; j < end; j++) {
    3006                                    int jRow = column[j];
    3007                                    int jNewRow = permuteInverse_[jRow];
    3008                                    if (jNewRow > iRow) {
    3009                                         work[jNewRow] = elementByRow[j];
    3010                                         largest = CoinMax(largest, CoinAbs(elementByRow[j]));
    3011                                    }
    3012                               }
    3013                               // slack - should it be permute
    3014                               kRow = permuteInverse_[kRow+numberColumns];
    3015                               if (kRow > iRow)
    3016                                    work[kRow] = -1.0;
    3017                          }
    3018                          CoinBigIndex j;
    3019                          int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    3020                          for (j = 0; j < number; j++) {
    3021                               int jRow = which[j];
    3022                               put[j] = work[jRow];
    3023                               work[jRow] = 0.0;
    3024                          }
    3025                     }
    3026                } else {
    3027                     // quadratic
    3028                     const int * columnQuadratic = quadratic->getIndices();
    3029                     const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    3030                     const int * columnQuadraticLength = quadratic->getVectorLengths();
    3031                     const double * quadraticElement = quadratic->getElements();
    3032                     for (iRow = 0; iRow < numberRows_; iRow++) {
    3033                          longDouble * put = sparseFactor_ + choleskyStart_[iRow];
    3034                          int * which = choleskyRow_ + indexStart_[iRow];
    3035                          int iOriginalRow = permute_[iRow];
    3036                          if (iOriginalRow < numberColumns) {
    3037                               CoinBigIndex j;
    3038                               iColumn = iOriginalRow;
    3039                               CoinWorkDouble value = diagonal[iColumn];
    3040                               if (CoinAbs(value) > 1.0e-100) {
    3041                                    value = 1.0 / value;
    3042                                    for (j = columnQuadraticStart[iColumn];
    3043                                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
    3044                                         int jColumn = columnQuadratic[j];
    3045                                         int jNewColumn = permuteInverse_[jColumn];
    3046                                         if (jNewColumn > iRow) {
    3047                                              work[jNewColumn] = -quadraticElement[j];
    3048                                         } else if (iColumn == jColumn) {
    3049                                              value += quadraticElement[j];
    3050                                         }
    3051                                    }
    3052                                    largest = CoinMax(largest, CoinAbs(value));
    3053                                    diagonal_[iRow] = -value;
    3054                                    CoinBigIndex start = columnStart[iColumn];
    3055                                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
    3056                                    for (j = start; j < end; j++) {
    3057                                         int kRow = row[j] + numberTotal;
    3058                                         kRow = permuteInverse_[kRow];
    3059                                         if (kRow > iRow) {
    3060                                              work[kRow] = element[j];
    3061                                              largest = CoinMax(largest, CoinAbs(element[j]));
    3062                                         }
    3063                                    }
    3064                               } else {
    3065                                    diagonal_[iRow] = -value;
    3066                               }
    3067                          } else if (iOriginalRow < numberTotal) {
    3068                               CoinWorkDouble value = diagonal[iOriginalRow];
    3069                               if (CoinAbs(value) > 1.0e-100) {
    3070                                    value = 1.0 / value;
    3071                                    largest = CoinMax(largest, CoinAbs(value));
    3072                               } else {
    3073                                    value = 1.0e100;
    3074                               }
    3075                               diagonal_[iRow] = -value;
    3076                               int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
    3077                               if (kRow > iRow)
    3078                                    work[kRow] = -1.0;
    3079                          } else {
    3080                               diagonal_[iRow] = delta2;
    3081                               int kRow = iOriginalRow - numberTotal;
    3082                               CoinBigIndex start = rowStart[kRow];
    3083                               CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
    3084                               for (CoinBigIndex j = start; j < end; j++) {
    3085                                    int jRow = column[j];
    3086                                    int jNewRow = permuteInverse_[jRow];
    3087                                    if (jNewRow > iRow) {
    3088                                         work[jNewRow] = elementByRow[j];
    3089                                         largest = CoinMax(largest, CoinAbs(elementByRow[j]));
    3090                                    }
    3091                               }
    3092                               // slack - should it be permute
    3093                               kRow = permuteInverse_[kRow+numberColumns];
    3094                               if (kRow > iRow)
    3095                                    work[kRow] = -1.0;
    3096                          }
    3097                          CoinBigIndex j;
    3098                          int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
    3099                          for (j = 0; j < number; j++) {
    3100                               int jRow = which[j];
    3101                               put[j] = work[jRow];
    3102                               work[jRow] = 0.0;
    3103                          }
    3104                          for (j = 0; j < numberRows_; j++)
    3105                               assert (!work[j]);
    3106                     }
    3107                }
     2770            denseDiagonal[iDense++] = 1.0;
     2771          }
     2772          dense += numberRows_;
     2773        }
     2774      }
     2775    }
     2776    CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
     2777    delta2 *= delta2;
     2778    // largest in initial matrix
     2779    CoinWorkDouble largest2 = 1.0e-20;
     2780    for (iRow = 0; iRow < numberRows_; iRow++) {
     2781      longDouble *put = sparseFactor_ + choleskyStart_[iRow];
     2782      int *which = choleskyRow_ + indexStart_[iRow];
     2783      int iOriginalRow = permute_[iRow];
     2784      int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
     2785      if (!rowLength[iOriginalRow])
     2786        rowsDropped_[iOriginalRow] = 1;
     2787      if (!rowsDropped_[iOriginalRow]) {
     2788        CoinBigIndex startRow = rowStart[iOriginalRow];
     2789        CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
     2790        work[iRow] = diagonalSlack[iOriginalRow] + delta2;
     2791        for (CoinBigIndex k = startRow; k < endRow; k++) {
     2792          int iColumn = column[k];
     2793          if (!whichDense_ || !whichDense_[iColumn]) {
     2794            CoinBigIndex start = columnStart[iColumn];
     2795            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
     2796            CoinWorkDouble multiplier = diagonal[iColumn] * elementByRow[k];
     2797            for (CoinBigIndex j = start; j < end; j++) {
     2798              int jRow = row[j];
     2799              int jNewRow = permuteInverse_[jRow];
     2800              if (jNewRow >= iRow && !rowsDropped_[jRow]) {
     2801                CoinWorkDouble value = element[j] * multiplier;
     2802                work[jNewRow] += value;
     2803              }
     2804            }
     2805          }