/* $Id: ClpCholeskyBase.cpp 1878 2012-08-30 15:43:19Z forrest $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

/*----------------------------------------------------------------------------*/
/*      Ordering code - courtesy of Anshul Gupta                              */
/*	(C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
/*----------------------------------------------------------------------------*/

/*  A compact no-frills Approximate Minimum Local Fill ordering code.

    References:

[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
    Edward Rothberg, SGI Manuscript, April 1996.
[2] An Approximate Minimum Degree Ordering Algorithm.
    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
    University of Florida, December 1994.
*/
/*----------------------------------------------------------------------------*/


#include "CoinPragma.hpp"

#include <iostream>

#include "ClpCholeskyBase.hpp"
#include "ClpInterior.hpp"
#include "ClpHelperFunctions.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinSort.hpp"
#include "ClpCholeskyDense.hpp"
#include "ClpMessage.hpp"
#include "ClpQuadraticObjective.hpp"

//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpCholeskyBase::ClpCholeskyBase (int denseThreshold) :
     type_(0),
     doKKT_(false),
     goDense_(0.7),
     choleskyCondition_(0.0),
     model_(NULL),
     numberTrials_(),
     numberRows_(0),
     status_(0),
     rowsDropped_(NULL),
     permuteInverse_(NULL),
     permute_(NULL),
     numberRowsDropped_(0),
     sparseFactor_(NULL),
     choleskyStart_(NULL),
     choleskyRow_(NULL),
     indexStart_(NULL),
     diagonal_(NULL),
     workDouble_(NULL),
     link_(NULL),
     workInteger_(NULL),
     clique_(NULL),
     sizeFactor_(0),
     sizeIndex_(0),
     firstDense_(0),
     rowCopy_(NULL),
     whichDense_(NULL),
     denseColumn_(NULL),
     dense_(NULL),
     denseThreshold_(denseThreshold)
{
     memset(integerParameters_, 0, 64 * sizeof(int));
     memset(doubleParameters_, 0, 64 * sizeof(double));
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpCholeskyBase::ClpCholeskyBase (const ClpCholeskyBase & rhs) :
     type_(rhs.type_),
     doKKT_(rhs.doKKT_),
     goDense_(rhs.goDense_),
     choleskyCondition_(rhs.choleskyCondition_),
     model_(rhs.model_),
     numberTrials_(rhs.numberTrials_),
     numberRows_(rhs.numberRows_),
     status_(rhs.status_),
     numberRowsDropped_(rhs.numberRowsDropped_)
{
     rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
     permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
     permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
     sizeFactor_ = rhs.sizeFactor_;
     sizeIndex_ = rhs.sizeIndex_;
     firstDense_ = rhs.firstDense_;
     sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
     choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
     indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
     choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
     diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
#if CLP_LONG_CHOLESKY!=1
     workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
#else
     // actually long double
     workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_), numberRows_));
#endif
     link_ = ClpCopyOfArray(rhs.link_, numberRows_);
     workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
     clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
     CoinMemcpyN(rhs.integerParameters_, 64, integerParameters_);
     CoinMemcpyN(rhs.doubleParameters_, 64, doubleParameters_);
     rowCopy_ = rhs.rowCopy_->clone();
     whichDense_ = NULL;
     denseColumn_ = NULL;
     dense_ = NULL;
     denseThreshold_ = rhs.denseThreshold_;
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpCholeskyBase::~ClpCholeskyBase ()
{
     delete [] rowsDropped_;
     delete [] permuteInverse_;
     delete [] permute_;
     delete [] sparseFactor_;
     delete [] choleskyStart_;
     delete [] choleskyRow_;
     delete [] indexStart_;
     delete [] diagonal_;
     delete [] workDouble_;
     delete [] link_;
     delete [] workInteger_;
     delete [] clique_;
     delete rowCopy_;
     delete [] whichDense_;
     delete [] denseColumn_;
     delete dense_;
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpCholeskyBase &
ClpCholeskyBase::operator=(const ClpCholeskyBase& rhs)
{
     if (this != &rhs) {
          type_ = rhs.type_;
          doKKT_ = rhs.doKKT_;
          goDense_ = rhs.goDense_;
          choleskyCondition_ = rhs.choleskyCondition_;
          model_ = rhs.model_;
          numberTrials_ = rhs.numberTrials_;
          numberRows_ = rhs.numberRows_;
          status_ = rhs.status_;
          numberRowsDropped_ = rhs.numberRowsDropped_;
          delete [] rowsDropped_;
          delete [] permuteInverse_;
          delete [] permute_;
          delete [] sparseFactor_;
          delete [] choleskyStart_;
          delete [] choleskyRow_;
          delete [] indexStart_;
          delete [] diagonal_;
          delete [] workDouble_;
          delete [] link_;
          delete [] workInteger_;
          delete [] clique_;
          delete rowCopy_;
          delete [] whichDense_;
          delete [] denseColumn_;
          delete dense_;
          rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
          permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
          permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
          sizeFactor_ = rhs.sizeFactor_;
          sizeIndex_ = rhs.sizeIndex_;
          firstDense_ = rhs.firstDense_;
          sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
          choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
          choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, rhs.sizeFactor_);
          indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
          choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
          diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
#if CLP_LONG_CHOLESKY!=1
          workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
#else
          // actually long double
          workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_), numberRows_));
#endif
          link_ = ClpCopyOfArray(rhs.link_, numberRows_);
          workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
          clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
          delete rowCopy_;
          rowCopy_ = rhs.rowCopy_->clone();
          whichDense_ = NULL;
          denseColumn_ = NULL;
          dense_ = NULL;
          denseThreshold_ = rhs.denseThreshold_;
     }
     return *this;
}
// reset numberRowsDropped and rowsDropped.
void
ClpCholeskyBase::resetRowsDropped()
{
     numberRowsDropped_ = 0;
     memset(rowsDropped_, 0, numberRows_);
}
/* Uses factorization to solve. - given as if KKT.
   region1 is rows+columns, region2 is rows */
void
ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
                           CoinWorkDouble diagonalScaleFactor)
{
     if (!doKKT_) {
          int iColumn;
          int numberColumns = model_->numberColumns();
          int numberTotal = numberRows_ + numberColumns;
          CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal];
          for (iColumn = 0; iColumn < numberTotal; iColumn++) {
               region1[iColumn] *= diagonal[iColumn];
               region1Save[iColumn] = region1[iColumn];
          }
          multiplyAdd(region1 + numberColumns, numberRows_, -1.0, region2, 1.0);
          model_->clpMatrix()->times(1.0, region1, region2);
          CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
          CoinWorkDouble scale = 1.0;
          CoinWorkDouble unscale = 1.0;
          if (maximumRHS > 1.0e-30) {
               if (maximumRHS <= 0.5) {
                    CoinWorkDouble factor = 2.0;
                    while (maximumRHS <= 0.5) {
                         maximumRHS *= factor;
                         scale *= factor;
                    } /* endwhile */
               } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
                    CoinWorkDouble factor = 0.5;
                    while (maximumRHS >= 2.0) {
                         maximumRHS *= factor;
                         scale *= factor;
                    } /* endwhile */
               }
               unscale = diagonalScaleFactor / scale;
          } else {
               //effectively zero
               scale = 0.0;
               unscale = 0.0;
          }
          multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
          solve(region2);
          multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
          multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns, 0.0);
          CoinZeroN(region1, numberColumns);
          model_->clpMatrix()->transposeTimes(1.0, region2, region1);
          for (iColumn = 0; iColumn < numberTotal; iColumn++)
               region1[iColumn] = region1[iColumn] * diagonal[iColumn] - region1Save[iColumn];
          delete [] region1Save;
     } else {
          // KKT
          int numberRowsModel = model_->numberRows();
          int numberColumns = model_->numberColumns();
          int numberTotal = numberColumns + numberRowsModel;
          CoinWorkDouble * array = new CoinWorkDouble [numberRows_];
          CoinMemcpyN(region1, numberTotal, array);
          CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
          assert (numberRows_ >= numberRowsModel + numberTotal);
          solve(array);
          int iRow;
          for (iRow = 0; iRow < numberTotal; iRow++) {
               if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
		 COIN_DETAIL_PRINT(printf("row region1 %d dropped %g\n", iRow, array[iRow]));
               }
          }
          for (; iRow < numberRows_; iRow++) {
               if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
		 COIN_DETAIL_PRINT(printf("row region2 %d dropped %g\n", iRow, array[iRow]));
               }
          }
          CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
          CoinMemcpyN(array, numberTotal, region1);
          delete [] array;
     }
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpCholeskyBase * ClpCholeskyBase::clone() const
{
     return new ClpCholeskyBase(*this);
}
// Forms ADAT - returns nonzero if not enough memory
int
ClpCholeskyBase::preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT)
{
     delete rowCopy_;
     rowCopy_ = model_->clpMatrix()->reverseOrderedCopy();
     if (!doKKT) {
          numberRows_ = model_->numberRows();
          rowsDropped_ = new char [numberRows_];
          memset(rowsDropped_, 0, numberRows_);
          numberRowsDropped_ = 0;
          // Space for starts
          choleskyStart_ = new CoinBigIndex[numberRows_+1];
          const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
          const int * columnLength = model_->clpMatrix()->getVectorLengths();
          const int * row = model_->clpMatrix()->getIndices();
          const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
          const int * rowLength = rowCopy_->getVectorLengths();
          const int * column = rowCopy_->getIndices();
          // We need two arrays for counts
          int * which = new int [numberRows_];
          int * used = new int[numberRows_+1];
          CoinZeroN(used, numberRows_);
          int iRow;
          sizeFactor_ = 0;
          int numberColumns = model_->numberColumns();
          int numberDense = 0;
          //denseThreshold_=3;
          if (denseThreshold_ > 0) {
               delete [] whichDense_;
               delete [] denseColumn_;
               delete dense_;
               whichDense_ = new char[numberColumns];
               int iColumn;
               used[numberRows_] = 0;
               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                    int length = columnLength[iColumn];
                    used[length] += 1;
               }
               int nLong = 0;
               int stop = CoinMax(denseThreshold_ / 2, 100);
               for (iRow = numberRows_; iRow >= stop; iRow--) {
                    if (used[iRow])
		      COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
                    nLong += used[iRow];
                    if (nLong > 50 || nLong > (numberColumns >> 2))
                         break;
               }
               CoinZeroN(used, numberRows_);
               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                    if (columnLength[iColumn] < denseThreshold_) {
                         whichDense_[iColumn] = 0;
                    } else {
                         whichDense_[iColumn] = 1;
                         numberDense++;
                    }
               }
               if (!numberDense || numberDense > 100) {
                    // free
                    delete [] whichDense_;
                    whichDense_ = NULL;
                    denseColumn_ = NULL;
                    dense_ = NULL;
               } else {
                    // space for dense columns
                    denseColumn_ = new longDouble [numberDense*numberRows_];
                    // dense cholesky
                    dense_ = new ClpCholeskyDense();
                    dense_->reserveSpace(NULL, numberDense);
                    COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
               }
          }
          int offset = includeDiagonal ? 0 : 1;
          if (lowerTriangular)
               offset = -offset;
          for (iRow = 0; iRow < numberRows_; iRow++) {
               int number = 0;
               // make sure diagonal exists if includeDiagonal
               if (!offset) {
                    which[0] = iRow;
                    used[iRow] = 1;
                    number = 1;
               }
               CoinBigIndex startRow = rowStart[iRow];
               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
               if (lowerTriangular) {
                    for (CoinBigIndex k = startRow; k < endRow; k++) {
                         int iColumn = column[k];
                         if (!whichDense_ || !whichDense_[iColumn]) {
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = row[j];
                                   if (jRow <= iRow + offset) {
                                        if (!used[jRow]) {
                                             used[jRow] = 1;
                                             which[number++] = jRow;
                                        }
                                   }
                              }
                         }
                    }
               } else {
                    for (CoinBigIndex k = startRow; k < endRow; k++) {
                         int iColumn = column[k];
                         if (!whichDense_ || !whichDense_[iColumn]) {
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = row[j];
                                   if (jRow >= iRow + offset) {
                                        if (!used[jRow]) {
                                             used[jRow] = 1;
                                             which[number++] = jRow;
                                        }
                                   }
                              }
                         }
                    }
               }
               sizeFactor_ += number;
               int j;
               for (j = 0; j < number; j++)
                    used[which[j]] = 0;
          }
          delete [] which;
          // Now we have size - create arrays and fill in
          try {
               choleskyRow_ = new int [sizeFactor_];
          } catch (...) {
               // no memory
               delete [] choleskyStart_;
               choleskyStart_ = NULL;
               return -1;
          }
          sizeFactor_ = 0;
          which = choleskyRow_;
          for (iRow = 0; iRow < numberRows_; iRow++) {
               int number = 0;
               // make sure diagonal exists if includeDiagonal
               if (!offset) {
                    which[0] = iRow;
                    used[iRow] = 1;
                    number = 1;
               }
               choleskyStart_[iRow] = sizeFactor_;
               CoinBigIndex startRow = rowStart[iRow];
               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
               if (lowerTriangular) {
                    for (CoinBigIndex k = startRow; k < endRow; k++) {
                         int iColumn = column[k];
                         if (!whichDense_ || !whichDense_[iColumn]) {
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = row[j];
                                   if (jRow <= iRow + offset) {
                                        if (!used[jRow]) {
                                             used[jRow] = 1;
                                             which[number++] = jRow;
                                        }
                                   }
                              }
                         }
                    }
               } else {
                    for (CoinBigIndex k = startRow; k < endRow; k++) {
                         int iColumn = column[k];
                         if (!whichDense_ || !whichDense_[iColumn]) {
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = row[j];
                                   if (jRow >= iRow + offset) {
                                        if (!used[jRow]) {
                                             used[jRow] = 1;
                                             which[number++] = jRow;
                                        }
                                   }
                              }
                         }
                    }
               }
               sizeFactor_ += number;
               int j;
               for (j = 0; j < number; j++)
                    used[which[j]] = 0;
               // Sort
               std::sort(which, which + number);
               // move which on
               which += number;
          }
          choleskyStart_[numberRows_] = sizeFactor_;
          delete [] used;
          return 0;
     } else {
          int numberRowsModel = model_->numberRows();
          int numberColumns = model_->numberColumns();
          int numberTotal = numberColumns + numberRowsModel;
          numberRows_ = 2 * numberRowsModel + numberColumns;
          rowsDropped_ = new char [numberRows_];
          memset(rowsDropped_, 0, numberRows_);
          numberRowsDropped_ = 0;
          CoinPackedMatrix * quadratic = NULL;
          ClpQuadraticObjective * quadraticObj =
               (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
          if (quadraticObj)
               quadratic = quadraticObj->quadraticObjective();
          int numberElements = model_->clpMatrix()->getNumElements();
          numberElements = numberElements + 2 * numberRowsModel + numberTotal;
          if (quadratic)
               numberElements += quadratic->getNumElements();
          // Space for starts
          choleskyStart_ = new CoinBigIndex[numberRows_+1];
          const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
          const int * columnLength = model_->clpMatrix()->getVectorLengths();
          const int * row = model_->clpMatrix()->getIndices();
          //const double * element = model_->clpMatrix()->getElements();
          // Now we have size - create arrays and fill in
          try {
               choleskyRow_ = new int [numberElements];
          } catch (...) {
               // no memory
               delete [] choleskyStart_;
               choleskyStart_ = NULL;
               return -1;
          }
          int iRow, iColumn;

          sizeFactor_ = 0;
          // matrix
          if (lowerTriangular) {
               if (!quadratic) {
                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                         choleskyStart_[iColumn] = sizeFactor_;
                         choleskyRow_[sizeFactor_++] = iColumn;
                         CoinBigIndex start = columnStart[iColumn];
                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                         if (!includeDiagonal)
                              start++;
                         for (CoinBigIndex j = start; j < end; j++) {
                              choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
                         }
                    }
               } else {
                    // Quadratic
                    const int * columnQuadratic = quadratic->getIndices();
                    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
                    const int * columnQuadraticLength = quadratic->getVectorLengths();
                    //const double * quadraticElement = quadratic->getElements();
                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                         choleskyStart_[iColumn] = sizeFactor_;
                         if (includeDiagonal)
                              choleskyRow_[sizeFactor_++] = iColumn;
                         CoinBigIndex j;
                         for ( j = columnQuadraticStart[iColumn];
                                   j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
                              int jColumn = columnQuadratic[j];
                              if (jColumn > iColumn)
                                   choleskyRow_[sizeFactor_++] = jColumn;
                         }
                         CoinBigIndex start = columnStart[iColumn];
                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                         for ( j = start; j < end; j++) {
                              choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
                         }
                    }
               }
               // slacks
               for (; iColumn < numberTotal; iColumn++) {
                    choleskyStart_[iColumn] = sizeFactor_;
                    if (includeDiagonal)
                         choleskyRow_[sizeFactor_++] = iColumn;
                    choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
               }
               // Transpose - nonzero diagonal (may regularize)
               for (iRow = 0; iRow < numberRowsModel; iRow++) {
                    choleskyStart_[iRow+numberTotal] = sizeFactor_;
                    // diagonal
                    if (includeDiagonal)
                         choleskyRow_[sizeFactor_++] = iRow + numberTotal;
               }
               choleskyStart_[numberRows_] = sizeFactor_;
          } else {
               // transpose
               ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
               const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
               const int * rowLength = rowCopy->getVectorLengths();
               const int * column = rowCopy->getIndices();
               if (!quadratic) {
                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                         choleskyStart_[iColumn] = sizeFactor_;
                         if (includeDiagonal)
                              choleskyRow_[sizeFactor_++] = iColumn;
                    }
               } else {
                    // Quadratic
                    // transpose
                    CoinPackedMatrix quadraticT;
                    quadraticT.reverseOrderedCopyOf(*quadratic);
                    const int * columnQuadratic = quadraticT.getIndices();
                    const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
                    const int * columnQuadraticLength = quadraticT.getVectorLengths();
                    //const double * quadraticElement = quadraticT.getElements();
                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                         choleskyStart_[iColumn] = sizeFactor_;
                         for (CoinBigIndex j = columnQuadraticStart[iColumn];
                                   j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
                              int jColumn = columnQuadratic[j];
                              if (jColumn < iColumn)
                                   choleskyRow_[sizeFactor_++] = jColumn;
                         }
                         if (includeDiagonal)
                              choleskyRow_[sizeFactor_++] = iColumn;
                    }
               }
               int iRow;
               // slacks
               for (iRow = 0; iRow < numberRowsModel; iRow++) {
                    choleskyStart_[iRow+numberColumns] = sizeFactor_;
                    if (includeDiagonal)
                         choleskyRow_[sizeFactor_++] = iRow + numberColumns;
               }
               for (iRow = 0; iRow < numberRowsModel; iRow++) {
                    choleskyStart_[iRow+numberTotal] = sizeFactor_;
                    CoinBigIndex start = rowStart[iRow];
                    CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
                    for (CoinBigIndex j = start; j < end; j++) {
                         choleskyRow_[sizeFactor_++] = column[j];
                    }
                    // slack
                    choleskyRow_[sizeFactor_++] = numberColumns + iRow;
                    if (includeDiagonal)
                         choleskyRow_[sizeFactor_++] = iRow + numberTotal;
               }
               choleskyStart_[numberRows_] = sizeFactor_;
          }
     }
     return 0;
}
/* Orders rows and saves pointer to matrix.and model */
int
ClpCholeskyBase::order(ClpInterior * model)
{
     model_ = model;
#define BASE_ORDER 2
#if BASE_ORDER>0
     if (!doKKT_ && model_->numberRows() > 6) {
          if (preOrder(false, true, false))
               return -1;
          //rowsDropped_ = new char [numberRows_];
          numberRowsDropped_ = 0;
          memset(rowsDropped_, 0, numberRows_);
          //rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
          // approximate minimum degree
          return orderAMD();
     }
#endif
     int numberRowsModel = model_->numberRows();
     int numberColumns = model_->numberColumns();
     int numberTotal = numberColumns + numberRowsModel;
     CoinPackedMatrix * quadratic = NULL;
     ClpQuadraticObjective * quadraticObj =
          (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     if (quadraticObj)
          quadratic = quadraticObj->quadraticObjective();
     if (!doKKT_) {
          numberRows_ = model->numberRows();
     } else {
          numberRows_ = 2 * numberRowsModel + numberColumns;
     }
     rowsDropped_ = new char [numberRows_];
     numberRowsDropped_ = 0;
     memset(rowsDropped_, 0, numberRows_);
     rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     const int * columnLength = model_->clpMatrix()->getVectorLengths();
     const int * row = model_->clpMatrix()->getIndices();
     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
     const int * rowLength = rowCopy_->getVectorLengths();
     const int * column = rowCopy_->getIndices();
     // We need arrays for counts
     int * which = new int [numberRows_];
     int * used = new int[numberRows_+1];
     int * count = new int[numberRows_];
     CoinZeroN(count, numberRows_);
     CoinZeroN(used, numberRows_);
     int iRow;
     sizeFactor_ = 0;
     permute_ = new int[numberRows_];
     for (iRow = 0; iRow < numberRows_; iRow++)
          permute_[iRow] = iRow;
     if (!doKKT_) {
          int numberDense = 0;
          if (denseThreshold_ > 0) {
               delete [] whichDense_;
               delete [] denseColumn_;
               delete dense_;
               whichDense_ = new char[numberColumns];
               int iColumn;
               used[numberRows_] = 0;
               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                    int length = columnLength[iColumn];
                    used[length] += 1;
               }
               int nLong = 0;
               int stop = CoinMax(denseThreshold_ / 2, 100);
               for (iRow = numberRows_; iRow >= stop; iRow--) {
                    if (used[iRow])
		      COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
                    nLong += used[iRow];
                    if (nLong > 50 || nLong > (numberColumns >> 2))
                         break;
               }
               CoinZeroN(used, numberRows_);
               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                    if (columnLength[iColumn] < denseThreshold_) {
                         whichDense_[iColumn] = 0;
                    } else {
                         whichDense_[iColumn] = 1;
                         numberDense++;
                    }
               }
               if (!numberDense || numberDense > 100) {
                    // free
                    delete [] whichDense_;
                    whichDense_ = NULL;
                    denseColumn_ = NULL;
                    dense_ = NULL;
               } else {
                    // space for dense columns
                    denseColumn_ = new longDouble [numberDense*numberRows_];
                    // dense cholesky
                    dense_ = new ClpCholeskyDense();
                    dense_->reserveSpace(NULL, numberDense);
                    COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
               }
          }
          /*
             Get row counts and size
          */
          for (iRow = 0; iRow < numberRows_; iRow++) {
               int number = 1;
               // make sure diagonal exists
               which[0] = iRow;
               used[iRow] = 1;
               CoinBigIndex startRow = rowStart[iRow];
               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
               for (CoinBigIndex k = startRow; k < endRow; k++) {
                    int iColumn = column[k];
                    if (!whichDense_ || !whichDense_[iColumn]) {
                         CoinBigIndex start = columnStart[iColumn];
                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                         for (CoinBigIndex j = start; j < end; j++) {
                              int jRow = row[j];
                              if (jRow < iRow) {
                                   if (!used[jRow]) {
                                        used[jRow] = 1;
                                        which[number++] = jRow;
                                        count[jRow]++;
                                   }
                              }
                         }
                    }
               }
               sizeFactor_ += number;
               count[iRow] += number;
               int j;
               for (j = 0; j < number; j++)
                    used[which[j]] = 0;
          }
          CoinSort_2(count, count + numberRows_, permute_);
     } else {
          // KKT
          int numberElements = model_->clpMatrix()->getNumElements();
          numberElements = numberElements + 2 * numberRowsModel + numberTotal;
          if (quadratic)
               numberElements += quadratic->getNumElements();
          // off diagonal
          numberElements -= numberRows_;
          sizeFactor_ = numberElements;
          // If we sort we need to redo symbolic etc
     }
     delete [] which;
     delete [] used;
     delete [] count;
     permuteInverse_ = new int [numberRows_];
     for (iRow = 0; iRow < numberRows_; iRow++) {
          //permute_[iRow]=iRow; // force no permute
          //permute_[iRow]=numberRows_-1-iRow; // force odd permute
          //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
          permuteInverse_[permute_[iRow]] = iRow;
     }
     return 0;
}
#if BASE_ORDER==1
/* Orders rows and saves pointer to matrix.and model */
int
ClpCholeskyBase::orderAMD()
{
     permuteInverse_ = new int [numberRows_];
     permute_ = new int[numberRows_];
     // Do ordering
     int returnCode = 0;
     // get more space and full matrix
     int space = 6 * sizeFactor_ + 100000;
     int * temp = new int [space];
     int * which = new int[2*numberRows_];
     CoinBigIndex * tempStart = new CoinBigIndex [numberRows_+1];
     memset(which, 0, numberRows_ * sizeof(int));
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          which[iRow] += choleskyStart_[iRow+1] - choleskyStart_[iRow] - 1;
          for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
               int jRow = choleskyRow_[j];
               which[jRow]++;
          }
     }
     CoinBigIndex sizeFactor = 0;
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          int length = which[iRow];
          permute_[iRow] = length;
          tempStart[iRow] = sizeFactor;
          which[iRow] = sizeFactor;
          sizeFactor += length;
     }
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          assert (choleskyRow_[choleskyStart_[iRow]] == iRow);
          for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
               int jRow = choleskyRow_[j];
               int put = which[iRow];
               temp[put] = jRow;
               which[iRow]++;
               put = which[jRow];
               temp[put] = iRow;
               which[jRow]++;
          }
     }
     for (int iRow = 1; iRow < numberRows_; iRow++)
          assert (which[iRow-1] == tempStart[iRow]);
     CoinBigIndex lastSpace = sizeFactor;
     delete [] choleskyRow_;
     choleskyRow_ = temp;
     delete [] choleskyStart_;
     choleskyStart_ = tempStart;
     // linked lists of sizes and lengths
     int * first = new int [numberRows_];
     int * next = new int [numberRows_];
     int * previous = new int [numberRows_];
     char * mark = new char[numberRows_];
     memset(mark, 0, numberRows_);
     CoinBigIndex * sort = new CoinBigIndex [numberRows_];
     int * order = new int [numberRows_];
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          first[iRow] = -1;
          next[iRow] = -1;
          previous[iRow] = -1;
          permuteInverse_[iRow] = -1;
     }
     int large = 1000 + 2 * numberRows_;
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          int n = permute_[iRow];
          if (first[n] < 0) {
               first[n] = iRow;
               previous[iRow] = n + large;
               next[iRow] = n + 2 * large;
          } else {
               int k = first[n];
               assert (k < numberRows_);
               first[n] = iRow;
               previous[iRow] = n + large;
               assert (previous[k] == n + large);
               next[iRow] = k;
               previous[k] = iRow;
          }
     }
     int smallest = 0;
     int done = 0;
     int numberCompressions = 0;
     int numberExpansions = 0;
     while (smallest < numberRows_) {
          if (first[smallest] < 0 || first[smallest] > numberRows_) {
               smallest++;
               continue;
          }
          int iRow = first[smallest];
          if (false) {
               int kRow = -1;
               int ss = 999999;
               for (int jRow = numberRows_ - 1; jRow >= 0; jRow--) {
                    if (permuteInverse_[jRow] < 0) {
                         int length = permute_[jRow];
                         assert (length > 0);
                         for (CoinBigIndex j = choleskyStart_[jRow];
                                   j < choleskyStart_[jRow] + length; j++) {
                              int jjRow = choleskyRow_[j];
                              assert (permuteInverse_[jjRow] < 0);
                         }
                         if (length < ss) {
                              ss = length;
                              kRow = jRow;
                         }
                    }
               }
               assert (smallest == ss);
               printf("list chose %d - full chose %d - length %d\n",
                      iRow, kRow, ss);
          }
          int kNext = next[iRow];
          first[smallest] = kNext;
          if (kNext < numberRows_)
               previous[kNext] = previous[iRow];
          previous[iRow] = -1;
          next[iRow] = -1;
          permuteInverse_[iRow] = done;
          done++;
          // Now add edges
          CoinBigIndex start = choleskyStart_[iRow];
          CoinBigIndex end = choleskyStart_[iRow] + permute_[iRow];
          int nSave = 0;
          for (CoinBigIndex k = start; k < end; k++) {
               int kRow = choleskyRow_[k];
               assert (permuteInverse_[kRow] < 0);
               //if (permuteInverse_[kRow]<0)
               which[nSave++] = kRow;
          }
          for (int i = 0; i < nSave; i++) {
               int kRow = which[i];
               int length = permute_[kRow];
               CoinBigIndex start = choleskyStart_[kRow];
               CoinBigIndex end = choleskyStart_[kRow] + length;
               for (CoinBigIndex j = start; j < end; j++) {
                    int jRow = choleskyRow_[j];
                    mark[jRow] = 1;
               }
               mark[kRow] = 1;
               int n = nSave;
               for (int j = 0; j < nSave; j++) {
                    int jRow = which[j];
                    if (!mark[jRow]) {
                         which[n++] = jRow;
                    }
               }
               for (CoinBigIndex j = start; j < end; j++) {
                    int jRow = choleskyRow_[j];
                    mark[jRow] = 0;
               }
               mark[kRow] = 0;
               if (n > nSave) {
                    bool inPlace = (n - nSave == 1);
                    if (!inPlace) {
                         // extend
                         int length = n - nSave + end - start;
                         if (length + lastSpace > space) {
                              // need to compress
                              numberCompressions++;
                              int newN = 0;
                              for (int iRow = 0; iRow < numberRows_; iRow++) {
                                   CoinBigIndex start = choleskyStart_[iRow];
                                   if (permuteInverse_[iRow] < 0) {
                                        sort[newN] = start;
                                        order[newN++] = iRow;
                                   } else {
                                        choleskyStart_[iRow] = -1;
                                        permute_[iRow] = 0;
                                   }
                              }
                              CoinSort_2(sort, sort + newN, order);
                              sizeFactor = 0;
                              for (int k = 0; k < newN; k++) {
                                   int iRow = order[k];
                                   int length = permute_[iRow];
                                   CoinBigIndex start = choleskyStart_[iRow];
                                   choleskyStart_[iRow] = sizeFactor;
                                   for (int j = 0; j < length; j++)
                                        choleskyRow_[sizeFactor+j] = choleskyRow_[start+j];
                                   sizeFactor += length;
                              }
                              lastSpace = sizeFactor;
                              if ((sizeFactor + numberRows_ + 1000) * 4 > 3 * space) {
                                   // need to expand
                                   numberExpansions++;
                                   space = (3 * space) / 2;
                                   int * temp = new int [space];
                                   memcpy(temp, choleskyRow_, sizeFactor * sizeof(int));
                                   delete [] choleskyRow_;
                                   choleskyRow_ = temp;
                              }
                         }
                    }
                    // Now add
                    start = choleskyStart_[kRow];
                    end = choleskyStart_[kRow] + permute_[kRow];
                    if (!inPlace)
                         choleskyStart_[kRow] = lastSpace;
                    CoinBigIndex put = choleskyStart_[kRow];
                    for (CoinBigIndex j = start; j < end; j++) {
                         int jRow = choleskyRow_[j];
                         if (permuteInverse_[jRow] < 0)
                              choleskyRow_[put++] = jRow;
                    }
                    for (int j = nSave; j < n; j++) {
                         int jRow = which[j];
                         choleskyRow_[put++] = jRow;
                    }
                    if (!inPlace) {
                         permute_[kRow] = put - lastSpace;
                         lastSpace = put;
                    } else {
                         permute_[kRow] = put - choleskyStart_[kRow];
                    }
               } else {
                    // pack down for new counts
                    CoinBigIndex put = start;
                    for (CoinBigIndex j = start; j < end; j++) {
                         int jRow = choleskyRow_[j];
                         if (permuteInverse_[jRow] < 0)
                              choleskyRow_[put++] = jRow;
                    }
                    permute_[kRow] = put - start;
               }
               // take out
               int iNext = next[kRow];
               int iPrevious = previous[kRow];
               if (iPrevious < numberRows_) {
                    next[iPrevious] = iNext;
               } else {
                    assert (iPrevious == length + large);
                    first[length] = iNext;
               }
               if (iNext < numberRows_) {
                    previous[iNext] = iPrevious;
               } else {
                    assert (iNext == length + 2 * large);
               }
               // put in
               length = permute_[kRow];
               smallest = CoinMin(smallest, length);
               if (first[length] < 0 || first[length] > numberRows_) {
                    first[length] = kRow;
                    previous[kRow] = length + large;
                    next[kRow] = length + 2 * large;
               } else {
                    int k = first[length];
                    assert (k < numberRows_);
                    first[length] = kRow;
                    previous[kRow] = length + large;
                    assert (previous[k] == length + large);
                    next[kRow] = k;
                    previous[k] = kRow;
               }
          }
     }
     // do rest
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          if (permuteInverse_[iRow] < 0)
               permuteInverse_[iRow] = done++;
     }
     COIN_DETAIL_PRINT(printf("%d compressions, %d expansions\n",
			      numberCompressions, numberExpansions));
     assert (done == numberRows_);
     delete [] sort;
     delete [] order;
     delete [] which;
     delete [] mark;
     delete [] first;
     delete [] next;
     delete [] previous;
     delete [] choleskyRow_;
     choleskyRow_ = NULL;
     delete [] choleskyStart_;
     choleskyStart_ = NULL;
#ifndef NDEBUG
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          permute_[iRow] = -1;
     }
#endif
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          permute_[permuteInverse_[iRow]] = iRow;
     }
#ifndef NDEBUG
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          assert(permute_[iRow] >= 0 && permute_[iRow] < numberRows_);
     }
#endif
     return returnCode;
}
#elif BASE_ORDER==2
/*----------------------------------------------------------------------------*/
/*	(C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
/*----------------------------------------------------------------------------*/

/*  A compact no-frills Approximate Minimum Local Fill ordering code.

    References:

[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
    Edward Rothberg, SGI Manuscript, April 1996.
[2] An Approximate Minimum Degree Ordering Algorithm.
    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
    University of Florida, December 1994.
*/

#include <math.h>
#include <stdlib.h>

typedef int             WSI;

#define NORDTHRESH	7
#define MAXIW		2147000000

//#define WSSMP_DBG
#ifdef WSSMP_DBG
static void chk (WSI ind, WSI i, WSI lim)
{

     if (i <= 0 || i > lim) {
          printf("%d: chk: myamlf: WAR**: i, lim = %d, %d\n", ind, i, lim);
     }
}
#endif

static void
myamlf(WSI n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI varbl[],
       WSI snxt[], WSI perm[], WSI invp[], WSI head[], WSI lsize[],
       WSI flag[], WSI erscore[], WSI locaux, WSI adjln, WSI speed)
{
     WSI l, i, j, k;
     double x, y;
     WSI maxmum, fltag, nodeg, scln, nm1, deg, tn,
         locatns, ipp, jpp, nnode, numpiv, node,
         nodeln, currloc, counter, numii, mindeg,
         i0, i1, i2, i4, i5, i6, i7, i9,
         j0, j1, j2, j3, j4, j5, j6, j7, j8, j9;
     /*                                             n cannot be less than NORDTHRESH
      if (n < 3) {
         if (n > 0) perm[0] = invp[0] = 1;
         if (n > 1) perm[1] = invp[1] = 2;
         return;
      }
     */
#ifdef WSSMP_DBG
     printf("Myamlf: n, locaux, adjln, speed = %d,%d,%d,%d\n", n, locaux, adjln, speed);
     for (i = 0; i < n; i++) flag[i] = 0;
     k = xadj[0];
     for (i = 1; i <= n; i++) {
          j = k + dgree[i-1];
          if (j < k || k < 1) printf("ERR**: myamlf: i, j, k = %d,%d,%d\n", i, j, k);
          for (l = k; l < j; l++) {
               if (adjncy[l - 1] < 1 || adjncy[l - 1] > n || adjncy[l - 1] == i)
                    printf("ERR**: myamlf: i, l, adjj[l] = %d,%d,%d\n", i, l, adjncy[l - 1]);
               if (flag[adjncy[l - 1] - 1] == i)
                    printf("ERR**: myamlf: duplicate entry: %d - %d\n", i, adjncy[l - 1]);
               flag[adjncy[l - 1] - 1] = i;
               nm1 = adjncy[l - 1] - 1;
               for (fltag = xadj[nm1]; fltag < xadj[nm1] + dgree[nm1]; fltag++) {
                    if (adjncy[fltag - 1] == i) goto L100;
               }
               printf("ERR**: Unsymmetric %d %d\n", i, fltag);
L100:
               ;
          }
          k = j;
          if (k > locaux) printf("ERR**: i, j, k, locaux = %d, %d, %d, %d\n",
                                      i, j, k, locaux);
     }
     if (n*(n - 1) < locaux - 1) printf("ERR**: myamlf: too many nnz, n, ne = %d, %d\n",
                                             n, adjln);
     for (i = 0; i < n; i++) flag[i] = 1;
     if (n > 10000) printf ("Finished error checking in AMF\n");
     for (i = locaux; i <= adjln; i++) adjncy[i - 1] = -22;
#endif

     maxmum = 0;
     mindeg = 1;
     fltag = 2;
     locatns = locaux - 1;
     nm1 = n - 1;
     counter = 1;
     for (l = 0; l < n; l++) {
          j = erscore[l];
#ifdef WSSMP_DBG
          chk(1, j, n);
#endif
          if (j > 0) {
               nnode = head[j - 1];
               if (nnode) perm[nnode - 1] = l + 1;
               snxt[l] = nnode;
               head[j - 1] = l + 1;
          } else {
               invp[l] = -(counter++);
               flag[l] = xadj[l] = 0;
          }
     }
     while (counter <= n) {
          for (deg = mindeg; head[deg - 1] < 1; deg++) {};
          nodeg = 0;
#ifdef WSSMP_DBG
          chk(2, deg, n - 1);
#endif
          node = head[deg - 1];
#ifdef WSSMP_DBG
          chk(3, node, n);
#endif
          mindeg = deg;
          nnode = snxt[node - 1];
          if (nnode) perm[nnode - 1] = 0;
          head[deg - 1] = nnode;
          nodeln = invp[node - 1];
          numpiv = varbl[node - 1];
          invp[node - 1] = -counter;
          counter += numpiv;
          varbl[node - 1] = -numpiv;
          if (nodeln) {
               j4 = locaux;
               i5 = lsize[node - 1] - nodeln;
               i2 = nodeln + 1;
               l = xadj[node - 1];
               for (i6 = 1; i6 <= i2; ++i6) {
                    if (i6 == i2) {
                         tn = node, i0 = l, scln = i5;
                    } else {
#ifdef WSSMP_DBG
                         chk(4, l, adjln);
#endif
                         tn = adjncy[l-1];
                         l++;
#ifdef WSSMP_DBG
                         chk(5, tn, n);
#endif
                         i0 = xadj[tn - 1];
                         scln = lsize[tn - 1];
                    }
                    for (i7 = 1; i7 <= scln; ++i7) {
#ifdef WSSMP_DBG
                         chk(6, i0, adjln);
#endif
                         i = adjncy[i0 - 1];
                         i0++;
#ifdef WSSMP_DBG
                         chk(7, i, n);
#endif
                         numii = varbl[i - 1];
                         if (numii > 0) {
                              if (locaux > adjln) {
                                   lsize[node - 1] -= i6;
                                   xadj[node - 1] = l;
                                   if (!lsize[node - 1]) xadj[node - 1] = 0;
                                   xadj[tn - 1] = i0;
                                   lsize[tn - 1] = scln - i7;
                                   if (!lsize[tn - 1]) xadj[tn - 1] = 0;
                                   for (j = 1; j <= n; j++) {
                                        i4 = xadj[j - 1];
                                        if (i4 > 0) {
                                             xadj[j - 1] = adjncy[i4 - 1];
#ifdef WSSMP_DBG
                                             chk(8, i4, adjln);
#endif
                                             adjncy[i4 - 1] = -j;
                                        }
                                   }
                                   i9 = j4 - 1;
                                   j6 = j7 = 1;
                                   while (j6 <= i9) {
#ifdef WSSMP_DBG
                                        chk(9, j6, adjln);
#endif
                                        j = -adjncy[j6 - 1];
                                        j6++;
                                        if (j > 0) {
#ifdef WSSMP_DBG
                                             chk(10, j7, adjln);
                                             chk(11, j, n);
#endif
                                             adjncy[j7 - 1] = xadj[j - 1];
                                             xadj[j - 1] = j7++;
                                             j8 = lsize[j - 1] - 1 + j7;
                                             for (; j7 < j8; j7++, j6++) adjncy[j7-1] = adjncy[j6-1];
                                        }
                                   }
                                   for (j0 = j7; j4 < locaux; j4++, j7++) {
#ifdef WSSMP_DBG
                                        chk(12, j4, adjln);
#endif
                                        adjncy[j7 - 1] = adjncy[j4 - 1];
                                   }
                                   j4 = j0;
                                   locaux = j7;
                                   i0 = xadj[tn - 1];
                                   l = xadj[node - 1];
                              }
                              adjncy[locaux-1] = i;
                              locaux++;
                              varbl[i - 1] = -numii;
                              nodeg += numii;
                              ipp = perm[i - 1];
                              nnode = snxt[i - 1];
#ifdef WSSMP_DBG
                              if (ipp) chk(13, ipp, n);
                              if (nnode) chk(14, nnode, n);
#endif
                              if (ipp) snxt[ipp - 1] = nnode;
                              else head[erscore[i - 1] - 1] = nnode;
                              if (nnode) perm[nnode - 1] = ipp;
                         }
                    }
                    if (tn != node) {
                         flag[tn - 1] = 0, xadj[tn - 1] = -node;
                    }
               }
               currloc = (j5 = locaux) - j4;
               locatns += currloc;
          } else {
               i1 = (j4 = xadj[node - 1]) + lsize[node - 1];
               for (j = j5 = j4; j < i1; j++) {
#ifdef WSSMP_DBG
                    chk(15, j, adjln);
#endif
                    i = adjncy[j - 1];
#ifdef WSSMP_DBG
                    chk(16, i, n);
#endif
                    numii = varbl[i - 1];
                    if (numii > 0) {
                         nodeg += numii;
                         varbl[i - 1] = -numii;
#ifdef WSSMP_DBG
                         chk(17, j5, adjln);
#endif
                         adjncy[j5-1] = i;
                         ipp = perm[i - 1];
                         nnode = snxt[i - 1];
                         j5++;
#ifdef WSSMP_DBG
                         if (ipp) chk(18, ipp, n);
                         if (nnode) chk(19, nnode, n);
#endif
                         if (ipp) snxt[ipp - 1] = nnode;
                         else head[erscore[i - 1] - 1] = nnode;
                         if (nnode) perm[nnode - 1] = ipp;
                    }
               }
               currloc = 0;
          }
#ifdef WSSMP_DBG
          chk(20, node, n);
#endif
          lsize[node - 1] = j5 - (xadj[node - 1] = j4);
          dgree[node - 1] = nodeg;
          if (fltag + n < 0 || fltag + n > MAXIW) {
               for (i = 1; i <= n; i++) if (flag[i - 1]) flag[i - 1] = 1;
               fltag = 2;
          }
          for (j3 = j4; j3 < j5; j3++) {
#ifdef WSSMP_DBG
               chk(21, j3, adjln);
#endif
               i = adjncy[j3 - 1];
#ifdef WSSMP_DBG
               chk(22, i, n);
#endif
               j = invp[i - 1];
               if (j > 0) {
                    i4 = fltag - (numii = -varbl[i - 1]);
                    i2 = xadj[i - 1] + j;
                    for (l = xadj[i - 1]; l < i2; l++) {
#ifdef WSSMP_DBG
                         chk(23, l, adjln);
#endif
                         tn = adjncy[l - 1];
#ifdef WSSMP_DBG
                         chk(24, tn, n);
#endif
                         j9 = flag[tn - 1];
                         if (j9 >= fltag) j9 -= numii;
                         else if (j9) j9 = dgree[tn - 1] + i4;
                         flag[tn - 1] = j9;
                    }
               }
          }
          for (j3 = j4; j3 < j5; j3++) {
#ifdef WSSMP_DBG
               chk(25, j3, adjln);
#endif
               i = adjncy[j3 - 1];
               i5 = deg = 0;
#ifdef WSSMP_DBG
               chk(26, i, n);
#endif
               j1 = (i4 = xadj[i - 1]) + invp[i - 1];
               for (l = j0 = i4; l < j1; l++) {
#ifdef WSSMP_DBG
                    chk(27, l, adjln);
#endif
                    tn = adjncy[l - 1];
#ifdef WSSMP_DBG
                    chk(70, tn, n);
#endif
                    j8 = flag[tn - 1];
                    if (j8) {
                         deg += (j8 - fltag);
#ifdef WSSMP_DBG
                         chk(28, i4, adjln);
#endif
                         adjncy[i4-1] = tn;
                         i5 += tn;
                         i4++;
                         while (i5 >= nm1) i5 -= nm1;
                    }
               }
               invp[i - 1] = (j2 = i4) - j0 + 1;
               i2 = j0 + lsize[i - 1];
               for (l = j1; l < i2; l++) {
#ifdef WSSMP_DBG
                    chk(29, l, adjln);
#endif
                    j = adjncy[l - 1];
#ifdef WSSMP_DBG
                    chk(30, j, n);
#endif
                    numii = varbl[j - 1];
                    if (numii > 0) {
                         deg += numii;
#ifdef WSSMP_DBG
                         chk(31, i4, adjln);
#endif
                         adjncy[i4-1] = j;
                         i5 += j;
                         i4++;
                         while (i5 >= nm1) i5 -= nm1;
                    }
               }
               if (invp[i - 1] == 1 && j2 == i4) {
                    numii = -varbl[i - 1];
                    xadj[i - 1] = -node;
                    nodeg -= numii;
                    counter += numii;
                    numpiv += numii;
                    invp[i - 1] = varbl[i - 1] = 0;
               } else {
                    if (dgree[i - 1] > deg) dgree[i - 1] = deg;
#ifdef WSSMP_DBG
                    chk(32, j2, adjln);
                    chk(33, j0, adjln);
#endif
                    adjncy[i4 - 1] = adjncy[j2 - 1];
                    adjncy[j2 - 1] = adjncy[j0 - 1];
                    adjncy[j0 - 1] = node;
                    lsize[i - 1] = i4 - j0 + 1;
                    i5++;
#ifdef WSSMP_DBG
                    chk(35, i5, n);
#endif
                    j = head[i5 - 1];
                    if (j > 0) {
                         snxt[i - 1] = perm[j - 1];
                         perm[j - 1] = i;
                    } else {
                         snxt[i - 1] = -j;
                         head[i5 - 1] = -i;
                    }
                    perm[i - 1] = i5;
               }
          }
#ifdef WSSMP_DBG
          chk(36, node, n);
#endif
          dgree[node - 1] = nodeg;
          if (maxmum < nodeg) maxmum = nodeg;
          fltag += maxmum;
#ifdef WSSMP_DBG
          if (fltag + n + 1 < 0) printf("ERR**: fltag = %d (A)\n", fltag);
#endif
          for (j3 = j4; j3 < j5; ++j3) {
#ifdef WSSMP_DBG
               chk(37, j3, adjln);
#endif
               i = adjncy[j3 - 1];
#ifdef WSSMP_DBG
               chk(38, i, n);
#endif
               if (varbl[i - 1] < 0) {
                    i5 = perm[i - 1];
#ifdef WSSMP_DBG
                    chk(39, i5, n);
#endif
                    j = head[i5 - 1];
                    if (j) {
                         if (j < 0) {
                              head[i5 - 1] = 0, i = -j;
                         } else {
#ifdef WSSMP_DBG
                              chk(40, j, n);
#endif
                              i = perm[j - 1];
                              perm[j - 1] = 0;
                         }
                         while (i) {
#ifdef WSSMP_DBG
                              chk(41, i, n);
#endif
                              if (!snxt[i - 1]) {
                                   i = 0;
                              } else {
                                   k = invp[i - 1];
                                   i2 = xadj[i - 1] + (scln = lsize[i - 1]);
                                   for (l = xadj[i - 1] + 1; l < i2; l++) {
#ifdef WSSMP_DBG
                                        chk(42, l, adjln);
                                        chk(43, adjncy[l - 1], n);
#endif
                                        flag[adjncy[l - 1] - 1] = fltag;
                                   }
                                   jpp = i;
                                   j = snxt[i - 1];
                                   while (j) {
#ifdef WSSMP_DBG
                                        chk(44, j, n);
#endif
                                        if (lsize[j - 1] == scln && invp[j - 1] == k) {
                                             i2 = xadj[j - 1] + scln;
                                             j8 = 1;
                                             for (l = xadj[j - 1] + 1; l < i2; l++) {
#ifdef WSSMP_DBG
                                                  chk(45, l, adjln);
                                                  chk(46, adjncy[l - 1], n);
#endif
                                                  if (flag[adjncy[l - 1] - 1] != fltag) {
                                                       j8 = 0;
                                                       break;
                                                  }
                                             }
                                             if (j8) {
                                                  xadj[j - 1] = -i;
                                                  varbl[i - 1] += varbl[j - 1];
                                                  varbl[j - 1] = invp[j - 1] = 0;
#ifdef WSSMP_DBG
                                                  chk(47, j, n);
                                                  chk(48, jpp, n);
#endif
                                                  j = snxt[j - 1];
                                                  snxt[jpp - 1] = j;
                                             } else {
                                                  jpp = j;
#ifdef WSSMP_DBG
                                                  chk(49, j, n);
#endif
                                                  j = snxt[j - 1];
                                             }
                                        } else {
                                             jpp = j;
#ifdef WSSMP_DBG
                                             chk(50, j, n);
#endif
                                             j = snxt[j - 1];
                                        }
                                   }
                                   fltag++;
#ifdef WSSMP_DBG
                                   if (fltag + n + 1 < 0) printf("ERR**: fltag = %d (B)\n", fltag);
#endif
#ifdef WSSMP_DBG
                                   chk(51, i, n);
#endif
                                   i = snxt[i - 1];
                              }
                         }
                    }
               }
          }
          j8 = nm1 - counter;
          switch (speed) {
          case 1:
               for (j = j3 = j4; j3 < j5; j3++) {
#ifdef WSSMP_DBG
                    chk(52, j3, adjln);
#endif
                    i = adjncy[j3 - 1];
#ifdef WSSMP_DBG
                    chk(53, i, n);
#endif
                    numii = varbl[i - 1];
                    if (numii < 0) {
                         k = 0;
                         i4 = (l = xadj[i - 1]) + invp[i - 1];
                         for (; l < i4; l++) {
#ifdef WSSMP_DBG
                              chk(54, l, adjln);
                              chk(55, adjncy[l - 1], n);
#endif
                              i5 = dgree[adjncy[l - 1] - 1];
                              if (k < i5) k = i5;
                         }
                         x = static_cast<double> (k - 1);
                         varbl[i - 1] = abs(numii);
                         j9 = dgree[i - 1] + nodeg;
                         deg = (j8 > j9 ? j9 : j8) + numii;
                         if (deg < 1) deg = 1;
                         y = static_cast<double> (deg);
                         dgree[i - 1] = deg;
                         /*
                                    printf("%i %f should >= %i %f\n",deg,y,k-1,x);
                                    if (y < x) printf("** \n"); else printf("\n");
                         */
                         y = y * y - y;
                         x = y - (x * x - x);
                         if (x < 1.1) x = 1.1;
                         deg = static_cast<WSI>(sqrt(x));
                         /*         if (deg < 1) deg = 1; */
                         erscore[i - 1] = deg;
#ifdef WSSMP_DBG
                         chk(56, deg, n - 1);
#endif
                         nnode = head[deg - 1];
                         if (nnode) perm[nnode - 1] = i;
                         snxt[i - 1] = nnode;
                         perm[i - 1] = 0;
#ifdef WSSMP_DBG
                         chk(57, j, adjln);
#endif
                         head[deg - 1] = adjncy[j-1] = i;
                         j++;
                         if (deg < mindeg) mindeg = deg;
                    }
               }
               break;
          case 2:
               for (j = j3 = j4; j3 < j5; j3++) {
#ifdef WSSMP_DBG
                    chk(58, j3, adjln);
#endif
                    i = adjncy[j3 - 1];
#ifdef WSSMP_DBG
                    chk(59, i, n);
#endif
                    numii = varbl[i - 1];
                    if (numii < 0) {
                         x = static_cast<double> (dgree[adjncy[xadj[i - 1] - 1] - 1] - 1);
                         varbl[i - 1] = abs(numii);
                         j9 = dgree[i - 1] + nodeg;
                         deg = (j8 > j9 ? j9 : j8) + numii;
                         if (deg < 1) deg = 1;
                         y = static_cast<double> (deg);
                         dgree[i - 1] = deg;
                         /*
                                    printf("%i %f should >= %i %f",deg,y,dgree[adjncy[xadj[i - 1] - 1] - 1]-1,x);
                                    if (y < x) printf("** \n"); else printf("\n");
                         */
                         y = y * y - y;
                         x = y - (x * x - x);
                         if (x < 1.1) x = 1.1;
                         deg = static_cast<WSI>(sqrt(x));
                         /*         if (deg < 1) deg = 1; */
                         erscore[i - 1] = deg;
#ifdef WSSMP_DBG
                         chk(60, deg, n - 1);
#endif
                         nnode = head[deg - 1];
                         if (nnode) perm[nnode - 1] = i;
                         snxt[i - 1] = nnode;
                         perm[i - 1] = 0;
#ifdef WSSMP_DBG
                         chk(61, j, adjln);
#endif
                         head[deg - 1] = adjncy[j-1] = i;
                         j++;
                         if (deg < mindeg) mindeg = deg;
                    }
               }
               break;
          default:
               for (j = j3 = j4; j3 < j5; j3++) {
#ifdef WSSMP_DBG
                    chk(62, j3, adjln);
#endif
                    i = adjncy[j3 - 1];
#ifdef WSSMP_DBG
                    chk(63, i, n);
#endif
                    numii = varbl[i - 1];
                    if (numii < 0) {
                         varbl[i - 1] = abs(numii);
                         j9 = dgree[i - 1] + nodeg;
                         deg = (j8 > j9 ? j9 : j8) + numii;
                         if (deg < 1) deg = 1;
                         dgree[i - 1] = deg;
#ifdef WSSMP_DBG
                         chk(64, deg, n - 1);
#endif
                         nnode = head[deg - 1];
                         if (nnode) perm[nnode - 1] = i;
                         snxt[i - 1] = nnode;
                         perm[i - 1] = 0;
#ifdef WSSMP_DBG
                         chk(65, j, adjln);
#endif
                         head[deg - 1] = adjncy[j-1] = i;
                         j++;
                         if (deg < mindeg) mindeg = deg;
                    }
               }
               break;
          }
          if (currloc) {
#ifdef WSSMP_DBG
               chk(66, node, n);
#endif
               locatns += (lsize[node - 1] - currloc), locaux = j;
          }
          varbl[node - 1] = numpiv + nodeg;
          lsize[node - 1] = j - j4;
          if (!lsize[node - 1]) flag[node - 1] = xadj[node - 1] = 0;
     }
     for (l = 1; l <= n; l++) {
          if (!invp[l - 1]) {
               for (i = -xadj[l - 1]; invp[i - 1] >= 0; i = -xadj[i - 1]) {};
               tn = i;
#ifdef WSSMP_DBG
               chk(67, tn, n);
#endif
               k = -invp[tn - 1];
               i = l;
#ifdef WSSMP_DBG
               chk(68, i, n);
#endif
               while (invp[i - 1] >= 0) {
                    nnode = -xadj[i - 1];
                    xadj[i - 1] = -tn;
                    if (!invp[i - 1]) invp[i - 1] = k++;
                    i = nnode;
               }
               invp[tn - 1] = -k;
          }
     }
     for (l = 0; l < n; l++) {
          i = abs(invp[l]);
#ifdef WSSMP_DBG
          chk(69, i, n);
#endif
          invp[l] = i;
          perm[i - 1] = l + 1;
     }
     return;
}
// This code is not needed, but left in in case needed sometime
#if 0
/*C--------------------------------------------------------------------------*/
void amlfdr(WSI *n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI *adjln,
            WSI *locaux, WSI varbl[], WSI snxt[], WSI perm[],
            WSI head[], WSI invp[], WSI lsize[], WSI flag[], WSI *ispeed)
{
     WSI nn, nlocaux, nadjln, speed, i, j, mx, mxj, *erscore;

#ifdef WSSMP_DBG
     printf("Calling amlfdr for n, speed = %d, %d\n", *n, *ispeed);
#endif

     if ((nn = *n) == 0) return;

#ifdef WSSMP_DBG
     if (nn == 31) {
          printf("n = %d; adjln = %d; locaux = %d; ispeed = %d\n",
                 *n, *adjln, *locaux, *ispeed);
     }
#endif

     if (nn < NORDTHRESH) {
          for (i = 0; i < nn; i++) lsize[i] = i;
          for (i = nn; i > 0; i--) {
               mx = dgree[0];
               mxj = 0;
               for (j = 1; j < i; j++)
                    if (dgree[j] > mx) {
                         mx = dgree[j];
                         mxj = j;
                    }
               invp[lsize[mxj]] = i;
               dgree[mxj] = dgree[i-1];
               lsize[mxj] = lsize[i-1];
          }
          return;
     }
     speed = *ispeed;
     if (speed < 3) {
          /*
              erscore = (WSI *)malloc(nn * sizeof(WSI));
              if (erscore == NULL) speed = 3;
          */
          wscmal (&nn, &i, &erscore);
          if (i != 0) speed = 3;
     }
     if (speed > 2) erscore = dgree;
     if (speed < 3) {
          for (i = 0; i < nn; i++) {
               perm[i] = 0;
               invp[i] = 0;
               head[i] = 0;
               flag[i] = 1;
               varbl[i] = 1;
               lsize[i] = dgree[i];
               erscore[i] = dgree[i];
          }
     } else {
          for (i = 0; i < nn; i++) {
               perm[i] = 0;
               invp[i] = 0;
               head[i] = 0;
               flag[i] = 1;
               varbl[i] = 1;
               lsize[i] = dgree[i];
          }
     }
     nlocaux = *locaux;
     nadjln = *adjln;

     myamlf(nn, xadj, adjncy, dgree, varbl, snxt, perm, invp,
            head, lsize, flag, erscore, nlocaux, nadjln, speed);
     /*
       if (speed < 3) free(erscore);
     */
     if (speed < 3) wscfr(&erscore);
     return;
}
#endif // end of taking out amlfdr
/*C--------------------------------------------------------------------------*/
#endif
// Orders rows
int
ClpCholeskyBase::orderAMD()
{
     permuteInverse_ = new int [numberRows_];
     permute_ = new int[numberRows_];
     // Do ordering
     int returnCode = 0;
     // get full matrix
     int space = 2 * sizeFactor_ + 10000 + 4 * numberRows_;
     int * temp = new int [space];
     CoinBigIndex * count = new CoinBigIndex [numberRows_];
     CoinBigIndex * tempStart = new CoinBigIndex [numberRows_+1];
     memset(count, 0, numberRows_ * sizeof(int));
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          count[iRow] += choleskyStart_[iRow+1] - choleskyStart_[iRow] - 1;
          for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
               int jRow = choleskyRow_[j];
               count[jRow]++;
          }
     }
#define OFFSET 1
     CoinBigIndex sizeFactor = 0;
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          int length = count[iRow];
          permute_[iRow] = length;
          // add 1 to starts
          tempStart[iRow] = sizeFactor + OFFSET;
          count[iRow] = sizeFactor;
          sizeFactor += length;
     }
     tempStart[numberRows_] = sizeFactor + OFFSET;
     // add 1 to rows
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          assert (choleskyRow_[choleskyStart_[iRow]] == iRow);
          for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow+1]; j++) {
               int jRow = choleskyRow_[j];
               int put = count[iRow];
               temp[put] = jRow + OFFSET;
               count[iRow]++;
               put = count[jRow];
               temp[put] = iRow + OFFSET;
               count[jRow]++;
          }
     }
     for (int iRow = 1; iRow < numberRows_; iRow++)
          assert (count[iRow-1] == tempStart[iRow] - OFFSET);
     delete [] choleskyRow_;
     choleskyRow_ = temp;
     delete [] choleskyStart_;
     choleskyStart_ = tempStart;
     int locaux = sizeFactor + OFFSET;
     delete [] count;
     int speed = integerParameters_[0];
     if (speed < 1 || speed > 2)
          speed = 3;
     int * use = new int [((speed<3) ? 7 : 6)*numberRows_];
     int * dgree = use;
     int * varbl = dgree + numberRows_;
     int * snxt = varbl + numberRows_;
     int * head = snxt + numberRows_;
     int * lsize = head + numberRows_;
     int * flag = lsize + numberRows_;
     int * erscore;
     for (int i = 0; i < numberRows_; i++) {
          dgree[i] = choleskyStart_[i+1] - choleskyStart_[i];
          head[i] = dgree[i];
          snxt[i] = 0;
          permute_[i] = 0;
          permuteInverse_[i] = 0;
          head[i] = 0;
          flag[i] = 1;
          varbl[i] = 1;
          lsize[i] = dgree[i];
     }
     if (speed < 3) {
          erscore = flag + numberRows_;
          for (int i = 0; i < numberRows_; i++)
               erscore[i] = dgree[i];
     } else {
          erscore = dgree;
     }
     myamlf(numberRows_, choleskyStart_, choleskyRow_,
            dgree, varbl, snxt, permute_, permuteInverse_,
            head, lsize, flag, erscore, locaux, space, speed);
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          permute_[iRow]--;
     }
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          permuteInverse_[permute_[iRow]] = iRow;
     }
     for (int iRow = 0; iRow < numberRows_; iRow++) {
          assert (permuteInverse_[iRow] >= 0 && permuteInverse_[iRow] < numberRows_);
     }
     delete [] use;
     delete [] choleskyRow_;
     choleskyRow_ = NULL;
     delete [] choleskyStart_;
     choleskyStart_ = NULL;
     return returnCode;
}
/* Does Symbolic factorization given permutation.
   This is called immediately after order.  If user provides this then
   user must provide factorize and solve.  Otherwise the default factorization is used
   returns non-zero if not enough memory */
int
ClpCholeskyBase::symbolic()
{
     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     const int * columnLength = model_->clpMatrix()->getVectorLengths();
     const int * row = model_->clpMatrix()->getIndices();
     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
     const int * rowLength = rowCopy_->getVectorLengths();
     const int * column = rowCopy_->getIndices();
     int numberRowsModel = model_->numberRows();
     int numberColumns = model_->numberColumns();
     int numberTotal = numberColumns + numberRowsModel;
     CoinPackedMatrix * quadratic = NULL;
     ClpQuadraticObjective * quadraticObj =
          (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     if (quadraticObj)
          quadratic = quadraticObj->quadraticObjective();
     // We need an array for counts
     int * used = new int[numberRows_+1];
     // If KKT then re-order so negative first
     if (doKKT_) {
          int nn = 0;
          int np = 0;
          int iRow;
          for (iRow = 0; iRow < numberRows_; iRow++) {
               int originalRow = permute_[iRow];
               if (originalRow < numberTotal)
                    permute_[nn++] = originalRow;
               else
                    used[np++] = originalRow;
          }
          CoinMemcpyN(used, np, permute_ + nn);
          for (iRow = 0; iRow < numberRows_; iRow++)
               permuteInverse_[permute_[iRow]] = iRow;
     }
     CoinZeroN(used, numberRows_);
     int iRow;
     int iColumn;
     bool noMemory = false;
     CoinBigIndex * Astart = new CoinBigIndex[numberRows_+1];
     int * Arow = NULL;
     try {
          Arow = new int [sizeFactor_];
     } catch (...) {
          // no memory
          delete [] Astart;
          return -1;
     }
     choleskyStart_ = new int[numberRows_+1];
     link_ = new int[numberRows_];
     workInteger_ = new CoinBigIndex[numberRows_];
     indexStart_ = new CoinBigIndex[numberRows_];
     clique_ = new int[numberRows_];
     // Redo so permuted upper triangular
     sizeFactor_ = 0;
     int * which = Arow;
     if (!doKKT_) {
          for (iRow = 0; iRow < numberRows_; iRow++) {
               int number = 0;
               int iOriginalRow = permute_[iRow];
               Astart[iRow] = sizeFactor_;
               CoinBigIndex startRow = rowStart[iOriginalRow];
               CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
               for (CoinBigIndex k = startRow; k < endRow; k++) {
                    int iColumn = column[k];
                    if (!whichDense_ || !whichDense_[iColumn]) {
                         CoinBigIndex start = columnStart[iColumn];
                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                         for (CoinBigIndex j = start; j < end; j++) {
                              int jRow = row[j];
                              int jNewRow = permuteInverse_[jRow];
                              if (jNewRow < iRow) {
                                   if (!used[jNewRow]) {
                                        used[jNewRow] = 1;
                                        which[number++] = jNewRow;
                                   }
                              }
                         }
                    }
               }
               sizeFactor_ += number;
               int j;
               for (j = 0; j < number; j++)
                    used[which[j]] = 0;
               // Sort
               std::sort(which, which + number);
               // move which on
               which += number;
          }
     } else {
          // KKT
          // transpose
          ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
          const int * rowLength = rowCopy->getVectorLengths();
          const int * column = rowCopy->getIndices();
          // temp
          bool permuted = false;
          for (iRow = 0; iRow < numberRows_; iRow++) {
               if (permute_[iRow] != iRow) {
                    permuted = true;
                    break;
               }
          }
          if (permuted) {
               // Need to permute - ugly
               if (!quadratic) {
                    for (iRow = 0; iRow < numberRows_; iRow++) {
                         Astart[iRow] = sizeFactor_;
                         int iOriginalRow = permute_[iRow];
                         if (iOriginalRow < numberColumns) {
                              // A may be upper triangular by mistake
                              iColumn = iOriginalRow;
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int kRow = row[j] + numberTotal;
                                   kRow = permuteInverse_[kRow];
                                   if (kRow < iRow)
                                        Arow[sizeFactor_++] = kRow;
                              }
                         } else if (iOriginalRow < numberTotal) {
                              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
                              if (kRow < iRow)
                                   Arow[sizeFactor_++] = kRow;
                         } else {
                              int kRow = iOriginalRow - numberTotal;
                              CoinBigIndex start = rowStart[kRow];
                              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = column[j];
                                   int jNewRow = permuteInverse_[jRow];
                                   if (jNewRow < iRow)
                                        Arow[sizeFactor_++] = jNewRow;
                              }
                              // slack - should it be permute
                              kRow = permuteInverse_[kRow+numberColumns];
                              if (kRow < iRow)
                                   Arow[sizeFactor_++] = kRow;
                         }
                         // Sort
                         std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
                    }
               } else {
                    // quadratic
                    // transpose
                    CoinPackedMatrix quadraticT;
                    quadraticT.reverseOrderedCopyOf(*quadratic);
                    const int * columnQuadratic = quadraticT.getIndices();
                    const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
                    const int * columnQuadraticLength = quadraticT.getVectorLengths();
                    for (iRow = 0; iRow < numberRows_; iRow++) {
                         Astart[iRow] = sizeFactor_;
                         int iOriginalRow = permute_[iRow];
                         if (iOriginalRow < numberColumns) {
                              // Quadratic bit
                              CoinBigIndex j;
                              for ( j = columnQuadraticStart[iOriginalRow];
                                        j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
                                   int jColumn = columnQuadratic[j];
                                   int jNewColumn = permuteInverse_[jColumn];
                                   if (jNewColumn < iRow)
                                        Arow[sizeFactor_++] = jNewColumn;
                              }
                              // A may be upper triangular by mistake
                              iColumn = iOriginalRow;
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (j = start; j < end; j++) {
                                   int kRow = row[j] + numberTotal;
                                   kRow = permuteInverse_[kRow];
                                   if (kRow < iRow)
                                        Arow[sizeFactor_++] = kRow;
                              }
                         } else if (iOriginalRow < numberTotal) {
                              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
                              if (kRow < iRow)
                                   Arow[sizeFactor_++] = kRow;
                         } else {
                              int kRow = iOriginalRow - numberTotal;
                              CoinBigIndex start = rowStart[kRow];
                              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = column[j];
                                   int jNewRow = permuteInverse_[jRow];
                                   if (jNewRow < iRow)
                                        Arow[sizeFactor_++] = jNewRow;
                              }
                              // slack - should it be permute
                              kRow = permuteInverse_[kRow+numberColumns];
                              if (kRow < iRow)
                                   Arow[sizeFactor_++] = kRow;
                         }
                         // Sort
                         std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
                    }
               }
          } else {
               if (!quadratic) {
                    for (iRow = 0; iRow < numberRows_; iRow++) {
                         Astart[iRow] = sizeFactor_;
                    }
               } else {
                    // Quadratic
                    // transpose
                    CoinPackedMatrix quadraticT;
                    quadraticT.reverseOrderedCopyOf(*quadratic);
                    const int * columnQuadratic = quadraticT.getIndices();
                    const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
                    const int * columnQuadraticLength = quadraticT.getVectorLengths();
                    //const double * quadraticElement = quadraticT.getElements();
                    for (iRow = 0; iRow < numberColumns; iRow++) {
                         int iOriginalRow = permute_[iRow];
                         Astart[iRow] = sizeFactor_;
                         for (CoinBigIndex j = columnQuadraticStart[iOriginalRow];
                                   j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
                              int jColumn = columnQuadratic[j];
                              int jNewColumn = permuteInverse_[jColumn];
                              if (jNewColumn < iRow)
                                   Arow[sizeFactor_++] = jNewColumn;
                         }
                    }
               }
               int iRow;
               // slacks
               for (iRow = 0; iRow < numberRowsModel; iRow++) {
                    Astart[iRow+numberColumns] = sizeFactor_;
               }
               for (iRow = 0; iRow < numberRowsModel; iRow++) {
                    Astart[iRow+numberTotal] = sizeFactor_;
                    CoinBigIndex start = rowStart[iRow];
                    CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
                    for (CoinBigIndex j = start; j < end; j++) {
                         Arow[sizeFactor_++] = column[j];
                    }
                    // slack
                    Arow[sizeFactor_++] = numberColumns + iRow;
               }
          }
          delete rowCopy;
     }
     Astart[numberRows_] = sizeFactor_;
     firstDense_ = numberRows_;
     symbolic1(Astart, Arow);
     // Now fill in indices
     try {
          // too big
          choleskyRow_ = new int[sizeFactor_];
     } catch (...) {
          // no memory
          noMemory = true;
     }
     double sizeFactor = sizeFactor_;
     if (!noMemory) {
          // Do lower triangular
          sizeFactor_ = 0;
          int * which = Arow;
          if (!doKKT_) {
               for (iRow = 0; iRow < numberRows_; iRow++) {
                    int number = 0;
                    int iOriginalRow = permute_[iRow];
                    Astart[iRow] = sizeFactor_;
                    if (!rowsDropped_[iOriginalRow]) {
                         CoinBigIndex startRow = rowStart[iOriginalRow];
                         CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
                         for (CoinBigIndex k = startRow; k < endRow; k++) {
                              int iColumn = column[k];
                              if (!whichDense_ || !whichDense_[iColumn]) {
                                   CoinBigIndex start = columnStart[iColumn];
                                   CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                                   for (CoinBigIndex j = start; j < end; j++) {
                                        int jRow = row[j];
                                        int jNewRow = permuteInverse_[jRow];
                                        if (jNewRow > iRow && !rowsDropped_[jRow]) {
                                             if (!used[jNewRow]) {
                                                  used[jNewRow] = 1;
                                                  which[number++] = jNewRow;
                                             }
                                        }
                                   }
                              }
                         }
                         sizeFactor_ += number;
                         int j;
                         for (j = 0; j < number; j++)
                              used[which[j]] = 0;
                         // Sort
                         std::sort(which, which + number);
                         // move which on
                         which += number;
                    }
               }
          } else {
               // KKT
               // temp
               bool permuted = false;
               for (iRow = 0; iRow < numberRows_; iRow++) {
                    if (permute_[iRow] != iRow) {
                         permuted = true;
                         break;
                    }
               }
               // but fake it
               for (iRow = 0; iRow < numberRows_; iRow++) {
                    //permute_[iRow]=iRow; // force no permute
                    //permuteInverse_[permute_[iRow]]=iRow;
               }
               if (permuted) {
                    // Need to permute - ugly
                    if (!quadratic) {
                         for (iRow = 0; iRow < numberRows_; iRow++) {
                              Astart[iRow] = sizeFactor_;
                              int iOriginalRow = permute_[iRow];
                              if (iOriginalRow < numberColumns) {
                                   iColumn = iOriginalRow;
                                   CoinBigIndex start = columnStart[iColumn];
                                   CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                                   for (CoinBigIndex j = start; j < end; j++) {
                                        int kRow = row[j] + numberTotal;
                                        kRow = permuteInverse_[kRow];
                                        if (kRow > iRow)
                                             Arow[sizeFactor_++] = kRow;
                                   }
                              } else if (iOriginalRow < numberTotal) {
                                   int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
                                   if (kRow > iRow)
                                        Arow[sizeFactor_++] = kRow;
                              } else {
                                   int kRow = iOriginalRow - numberTotal;
                                   CoinBigIndex start = rowStart[kRow];
                                   CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
                                   for (CoinBigIndex j = start; j < end; j++) {
                                        int jRow = column[j];
                                        int jNewRow = permuteInverse_[jRow];
                                        if (jNewRow > iRow)
                                             Arow[sizeFactor_++] = jNewRow;
                                   }
                                   // slack - should it be permute
                                   kRow = permuteInverse_[kRow+numberColumns];
                                   if (kRow > iRow)
                                        Arow[sizeFactor_++] = kRow;
                              }
                              // Sort
                              std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
                         }
                    } else {
                         // quadratic
                         const int * columnQuadratic = quadratic->getIndices();
                         const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
                         const int * columnQuadraticLength = quadratic->getVectorLengths();
                         for (iRow = 0; iRow < numberRows_; iRow++) {
                              Astart[iRow] = sizeFactor_;
                              int iOriginalRow = permute_[iRow];
                              if (iOriginalRow < numberColumns) {
                                   // Quadratic bit
                                   CoinBigIndex j;
                                   for ( j = columnQuadraticStart[iOriginalRow];
                                             j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
                                        int jColumn = columnQuadratic[j];
                                        int jNewColumn = permuteInverse_[jColumn];
                                        if (jNewColumn > iRow)
                                             Arow[sizeFactor_++] = jNewColumn;
                                   }
                                   iColumn = iOriginalRow;
                                   CoinBigIndex start = columnStart[iColumn];
                                   CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                                   for (j = start; j < end; j++) {
                                        int kRow = row[j] + numberTotal;
                                        kRow = permuteInverse_[kRow];
                                        if (kRow > iRow)
                                             Arow[sizeFactor_++] = kRow;
                                   }
                              } else if (iOriginalRow < numberTotal) {
                                   int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
                                   if (kRow > iRow)
                                        Arow[sizeFactor_++] = kRow;
                              } else {
                                   int kRow = iOriginalRow - numberTotal;
                                   CoinBigIndex start = rowStart[kRow];
                                   CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
                                   for (CoinBigIndex j = start; j < end; j++) {
                                        int jRow = column[j];
                                        int jNewRow = permuteInverse_[jRow];
                                        if (jNewRow > iRow)
                                             Arow[sizeFactor_++] = jNewRow;
                                   }
                                   // slack - should it be permute
                                   kRow = permuteInverse_[kRow+numberColumns];
                                   if (kRow > iRow)
                                        Arow[sizeFactor_++] = kRow;
                              }
                              // Sort
                              std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
                         }
                    }
               } else {
                    if (!quadratic) {
                         for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                              Astart[iColumn] = sizeFactor_;
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   Arow[sizeFactor_++] = row[j] + numberTotal;
                              }
                         }
                    } else {
                         // Quadratic
                         const int * columnQuadratic = quadratic->getIndices();
                         const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
                         const int * columnQuadraticLength = quadratic->getVectorLengths();
                         //const double * quadraticElement = quadratic->getElements();
                         for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                              Astart[iColumn] = sizeFactor_;
                              CoinBigIndex j;
                              for ( j = columnQuadraticStart[iColumn];
                                        j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
                                   int jColumn = columnQuadratic[j];
                                   if (jColumn > iColumn)
                                        Arow[sizeFactor_++] = jColumn;
                              }
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for ( j = start; j < end; j++) {
                                   Arow[sizeFactor_++] = row[j] + numberTotal;
                              }
                         }
                    }
                    // slacks
                    for (iRow = 0; iRow < numberRowsModel; iRow++) {
                         Astart[iRow+numberColumns] = sizeFactor_;
                         Arow[sizeFactor_++] = iRow + numberTotal;
                    }
                    // Transpose - nonzero diagonal (may regularize)
                    for (iRow = 0; iRow < numberRowsModel; iRow++) {
                         Astart[iRow+numberTotal] = sizeFactor_;
                    }
               }
               Astart[numberRows_] = sizeFactor_;
          }
          symbolic2(Astart, Arow);
          if (sizeIndex_ < sizeFactor_) {
               int * indices = NULL;
               try {
                    indices = new int[sizeIndex_];
               } catch (...) {
                    // no memory
                    noMemory = true;
               }
               if (!noMemory)  {
                    CoinMemcpyN(choleskyRow_, sizeIndex_, indices);
                    delete [] choleskyRow_;
                    choleskyRow_ = indices;
               }
          }
     }
     delete [] used;
     // Use cholesky regions
     delete [] Astart;
     delete [] Arow;
     double flops = 0.0;
     for (iRow = 0; iRow < numberRows_; iRow++) {
          int length = choleskyStart_[iRow+1] - choleskyStart_[iRow];
          flops += static_cast<double> (length) * (length + 2.0);
     }
     if (model_->messageHandler()->logLevel() > 0)
          std::cout << sizeFactor << " elements in sparse Cholesky, flop count " << flops << std::endl;
     try {
          sparseFactor_ = new longDouble [sizeFactor_];
#if CLP_LONG_CHOLESKY!=1
          workDouble_ = new longDouble[numberRows_];
#else
          // actually long double
          workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]);
#endif
          diagonal_ = new longDouble[numberRows_];
     } catch (...) {
          // no memory
          noMemory = true;
     }
     if (noMemory) {
          delete [] choleskyRow_;
          choleskyRow_ = NULL;
          delete [] choleskyStart_;
          choleskyStart_ = NULL;
          delete [] permuteInverse_;
          permuteInverse_ = NULL;
          delete [] permute_;
          permute_ = NULL;
          delete [] choleskyStart_;
          choleskyStart_ = NULL;
          delete [] indexStart_;
          indexStart_ = NULL;
          delete [] link_;
          link_ = NULL;
          delete [] workInteger_;
          workInteger_ = NULL;
          delete [] sparseFactor_;
          sparseFactor_ = NULL;
          delete [] workDouble_;
          workDouble_ = NULL;
          delete [] diagonal_;
          diagonal_ = NULL;
          delete [] clique_;
          clique_ = NULL;
          return -1;
     }
     return 0;
}
int
ClpCholeskyBase::symbolic1(const CoinBigIndex * Astart, const int * Arow)
{
     int * marked = reinterpret_cast<int *> (workInteger_);
     int iRow;
     // may not need to do this here but makes debugging easier
     for (iRow = 0; iRow < numberRows_; iRow++) {
          marked[iRow] = -1;
          link_[iRow] = -1;
          choleskyStart_[iRow] = 0; // counts
     }
     for (iRow = 0; iRow < numberRows_; iRow++) {
          marked[iRow] = iRow;
          for (CoinBigIndex j = Astart[iRow]; j < Astart[iRow+1]; j++) {
               int kRow = Arow[j];
               while (marked[kRow] != iRow ) {
                    if (link_[kRow] < 0 )
                         link_[kRow] = iRow;
                    choleskyStart_[kRow]++;
                    marked[kRow] = iRow;
                    kRow = link_[kRow];
               }
          }
     }
     sizeFactor_ = 0;
     for (iRow = 0; iRow < numberRows_; iRow++) {
          int number = choleskyStart_[iRow];
          choleskyStart_[iRow] = sizeFactor_;
          sizeFactor_ += number;
     }
     choleskyStart_[numberRows_] = sizeFactor_;
     return sizeFactor_;;
}
void
ClpCholeskyBase::symbolic2(const CoinBigIndex * Astart, const int * Arow)
{
     int * mergeLink = clique_;
     int * marker = reinterpret_cast<int *> (workInteger_);
     int iRow;
     for (iRow = 0; iRow < numberRows_; iRow++) {
          marker[iRow] = -1;
          mergeLink[iRow] = -1;
          link_[iRow] = -1; // not needed but makes debugging easier
     }
     int start = 0;
     int end = 0;
     choleskyStart_[0] = 0;

     for (iRow = 0; iRow < numberRows_; iRow++) {
          int nz = 0;
          int merge = mergeLink[iRow];
          bool marked = false;
          if (merge < 0)
               marker[iRow] = iRow;
          else
               marker[iRow] = merge;
          start = end;
          int startSub = start;
          link_[iRow] = numberRows_;
          CoinBigIndex j;
          for ( j = Astart[iRow]; j < Astart[iRow+1]; j++) {
               int kRow = Arow[j];
               int k = iRow;
               int linked = link_[iRow];
#ifndef NDEBUG
               int count = 0;
#endif
               while (linked <= kRow) {
                    k = linked;
                    linked = link_[k];
#ifndef NDEBUG
                    count++;
                    assert (count < numberRows_);
#endif
               }
               nz++;
               link_[k] = kRow;
               link_[kRow] = linked;
               if (marker[kRow] != marker[iRow])
                    marked = true;
          }
          bool reuse = false;
          // Check if we can re-use indices
          if (!marked && merge >= 0 && mergeLink[merge] < 0) {
               // can re-use all
               startSub = indexStart_[merge] + 1;
               nz = choleskyStart_[merge+1] - (choleskyStart_[merge] + 1);
               reuse = true;
          } else {
               // See if we can re-use any
               int k = mergeLink[iRow];
               int maxLength = 0;
               while (k >= 0) {
                    int length = choleskyStart_[k+1] - (choleskyStart_[k] + 1);
                    int start = indexStart_[k] + 1;
                    int stop = start + length;
                    if (length > maxLength) {
                         maxLength = length;
                         startSub = start;
                    }
                    int linked = iRow;
                    for (CoinBigIndex j = start; j < stop; j++) {
                         int kRow = choleskyRow_[j];
                         int kk = linked;
                         linked = link_[kk];
                         while (linked < kRow) {
                              kk = linked;
                              linked = link_[kk];
                         }
                         if (linked != kRow) {
                              nz++;
                              link_[kk] = kRow;
                              link_[kRow] = linked;
                              linked = kRow;
                         }
                    }
                    k = mergeLink[k];
               }
               if (nz == maxLength)
                    reuse = true; // can re-use
          }
          //reuse=false; //temp
          if (!reuse) {
               end += nz;
               startSub = start;
               int kRow = iRow;
               for (int j = start; j < end; j++) {
                    kRow = link_[kRow];
                    choleskyRow_[j] = kRow;
                    assert (kRow < numberRows_);
                    marker[kRow] = iRow;
               }
               marker[iRow] = iRow;
          }
          indexStart_[iRow] = startSub;
          choleskyStart_[iRow+1] = choleskyStart_[iRow] + nz;
          if (nz > 1) {
               int kRow = choleskyRow_[startSub];
               mergeLink[iRow] = mergeLink[kRow];
               mergeLink[kRow] = iRow;
          }
          // should not be needed
          //std::sort(choleskyRow_+indexStart_[iRow]
          //      ,choleskyRow_+indexStart_[iRow]+nz);
          //#define CLP_DEBUG
#ifdef CLP_DEBUG
          int last = -1;
          for ( j = indexStart_[iRow]; j < indexStart_[iRow] + nz; j++) {
               int kRow = choleskyRow_[j];
               assert (kRow > last);
               last = kRow;
          }
#endif
     }
     sizeFactor_ = choleskyStart_[numberRows_];
     sizeIndex_ = start;
     // find dense segment here
     int numberleft = numberRows_;
     for (iRow = 0; iRow < numberRows_; iRow++) {
          CoinBigIndex left = sizeFactor_ - choleskyStart_[iRow];
          double n = numberleft;
          double threshold = n * (n - 1.0) * 0.5 * goDense_;
          if ( left >= threshold)
               break;
          numberleft--;
     }
     //iRow=numberRows_;
     int nDense = numberRows_ - iRow;
#define DENSE_THRESHOLD 8
     // don't do if dense columns
     if (nDense >= DENSE_THRESHOLD && !dense_) {
       COIN_DETAIL_PRINT(printf("Going dense for last %d rows\n", nDense));
          // make sure we don't disturb any indices
          CoinBigIndex k = 0;
          for (int jRow = 0; jRow < iRow; jRow++) {
               int nz = choleskyStart_[jRow+1] - choleskyStart_[jRow];
               k = CoinMax(k, indexStart_[jRow] + nz);
          }
          indexStart_[iRow] = k;
          int j;
          for (j = iRow + 1; j < numberRows_; j++) {
               choleskyRow_[k++] = j;
               indexStart_[j] = k;
          }
          sizeIndex_ = k;
          assert (k <= sizeFactor_); // can't happen with any reasonable defaults
          k = choleskyStart_[iRow];
          for (j = iRow + 1; j <= numberRows_; j++) {
               k += numberRows_ - j;
               choleskyStart_[j] = k;
          }
          // allow for blocked dense
          ClpCholeskyDense dense;
          sizeFactor_ = choleskyStart_[iRow] + dense.space(nDense);
          firstDense_ = iRow;
          if (doKKT_) {
               // redo permute so negative ones first
               int putN = firstDense_;
               int putP = 0;
               int numberRowsModel = model_->numberRows();
               int numberColumns = model_->numberColumns();
               int numberTotal = numberColumns + numberRowsModel;
               for (iRow = firstDense_; iRow < numberRows_; iRow++) {
                    int originalRow = permute_[iRow];
                    if (originalRow < numberTotal)
                         permute_[putN++] = originalRow;
                    else
                         permuteInverse_[putP++] = originalRow;
               }
               for (iRow = putN; iRow < numberRows_; iRow++) {
                    permute_[iRow] = permuteInverse_[iRow-putN];
               }
               for (iRow = 0; iRow < numberRows_; iRow++) {
                    permuteInverse_[permute_[iRow]] = iRow;
               }
          }
     }
     // Clean up clique info
     for (iRow = 0; iRow < numberRows_; iRow++)
          clique_[iRow] = 0;
     int lastClique = -1;
     bool inClique = false;
     for (iRow = 1; iRow < firstDense_; iRow++) {
          int sizeLast = choleskyStart_[iRow] - choleskyStart_[iRow-1];
          int sizeThis = choleskyStart_[iRow+1] - choleskyStart_[iRow];
          if (indexStart_[iRow] == indexStart_[iRow-1] + 1 &&
                    sizeThis == sizeLast - 1 &&
                    sizeThis) {
               // in clique
               if (!inClique) {
                    inClique = true;
                    lastClique = iRow - 1;
               }
          } else if (inClique) {
               int sizeClique = iRow - lastClique;
               for (int i = lastClique; i < iRow; i++) {
                    clique_[i] = sizeClique;
                    sizeClique--;
               }
               inClique = false;
          }
     }
     if (inClique) {
          int sizeClique = iRow - lastClique;
          for (int i = lastClique; i < iRow; i++) {
               clique_[i] = sizeClique;
               sizeClique--;
          }
     }
     //for (iRow=0;iRow<numberRows_;iRow++)
     //clique_[iRow]=0;
}
/* Factorize - filling in rowsDropped and returning number dropped */
int
ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped)
{
     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     const int * columnLength = model_->clpMatrix()->getVectorLengths();
     const int * row = model_->clpMatrix()->getIndices();
     const double * element = model_->clpMatrix()->getElements();
     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
     const int * rowLength = rowCopy_->getVectorLengths();
     const int * column = rowCopy_->getIndices();
     const double * elementByRow = rowCopy_->getElements();
     int numberColumns = model_->clpMatrix()->getNumCols();
     //perturbation
     CoinWorkDouble perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
     //perturbation=perturbation*perturbation*100000000.0;
     if (perturbation > 1.0) {
#ifdef COIN_DEVELOP
          //if (model_->model()->logLevel()&4)
          std::cout << "large perturbation " << perturbation << std::endl;
#endif
          perturbation = CoinSqrt(perturbation);
          perturbation = 1.0;
     }
     int iRow;
     int iColumn;
     longDouble * work = workDouble_;
     CoinZeroN(work, numberRows_);
     int newDropped = 0;
     CoinWorkDouble largest = 1.0;
     CoinWorkDouble smallest = COIN_DBL_MAX;
     int numberDense = 0;
     if (!doKKT_) {
          const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
          if (dense_)
               numberDense = dense_->numberRows();
          if (whichDense_) {
               longDouble * denseDiagonal = dense_->diagonal();
               longDouble * dense = denseColumn_;
               int iDense = 0;
               for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                    if (whichDense_[iColumn]) {
                         CoinZeroN(dense, numberRows_);
                         CoinBigIndex start = columnStart[iColumn];
                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                         if (diagonal[iColumn]) {
                              denseDiagonal[iDense++] = 1.0 / diagonal[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = row[j];
                                   dense[jRow] = element[j];
                              }
                         } else {
                              denseDiagonal[iDense++] = 1.0;
                         }
                         dense += numberRows_;
                    }
               }
          }
          CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
          delta2 *= delta2;
          // largest in initial matrix
          CoinWorkDouble largest2 = 1.0e-20;
          for (iRow = 0; iRow < numberRows_; iRow++) {
               longDouble * put = sparseFactor_ + choleskyStart_[iRow];
               int * which = choleskyRow_ + indexStart_[iRow];
               int iOriginalRow = permute_[iRow];
               int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
               if (!rowLength[iOriginalRow])
                    rowsDropped_[iOriginalRow] = 1;
               if (!rowsDropped_[iOriginalRow]) {
                    CoinBigIndex startRow = rowStart[iOriginalRow];
                    CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
                    work[iRow] = diagonalSlack[iOriginalRow] + delta2;
                    for (CoinBigIndex k = startRow; k < endRow; k++) {
                         int iColumn = column[k];
                         if (!whichDense_ || !whichDense_[iColumn]) {
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              CoinWorkDouble multiplier = diagonal[iColumn] * elementByRow[k];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = row[j];
                                   int jNewRow = permuteInverse_[jRow];
                                   if (jNewRow >= iRow && !rowsDropped_[jRow]) {
                                        CoinWorkDouble value = element[j] * multiplier;
                                        work[jNewRow] += value;
                                   }
                              }
                         }
                    }
                    diagonal_[iRow] = work[iRow];
                    largest2 = CoinMax(largest2, CoinAbs(work[iRow]));
                    work[iRow] = 0.0;
                    int j;
                    for (j = 0; j < number; j++) {
                         int jRow = which[j];
                         put[j] = work[jRow];
                         largest2 = CoinMax(largest2, CoinAbs(work[jRow]));
                         work[jRow] = 0.0;
                    }
               } else {
                    // dropped
                    diagonal_[iRow] = 1.0;
                    int j;
                    for (j = 1; j < number; j++) {
                         put[j] = 0.0;
                    }
               }
          }
          //check sizes
          largest2 *= 1.0e-20;
          largest = CoinMin(largest2, CHOL_SMALL_VALUE);
          int numberDroppedBefore = 0;
          for (iRow = 0; iRow < numberRows_; iRow++) {
               int dropped = rowsDropped_[iRow];
               // Move to int array
               rowsDropped[iRow] = dropped;
               if (!dropped) {
                    CoinWorkDouble diagonal = diagonal_[iRow];
                    if (diagonal > largest2) {
                         diagonal_[iRow] = diagonal + perturbation;
                    } else {
                         diagonal_[iRow] = diagonal + perturbation;
                         rowsDropped[iRow] = 2;
                         numberDroppedBefore++;
                         //printf("dropped - small diagonal %g\n",diagonal);
                    }
               }
          }
          doubleParameters_[10] = CoinMax(1.0e-20, largest);
          integerParameters_[20] = 0;
          doubleParameters_[3] = 0.0;
          doubleParameters_[4] = COIN_DBL_MAX;
          integerParameters_[34] = 0; // say all must be positive
          factorizePart2(rowsDropped);
          newDropped = integerParameters_[20] + numberDroppedBefore;
          largest = doubleParameters_[3];
          smallest = doubleParameters_[4];
          if (model_->messageHandler()->logLevel() > 1)
               std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
          choleskyCondition_ = largest / smallest;
          if (whichDense_) {
               int i;
               for ( i = 0; i < numberRows_; i++) {
                    assert (diagonal_[i] >= 0.0);
                    diagonal_[i] = CoinSqrt(diagonal_[i]);
               }
               // Update dense columns (just L)
               // Zero out dropped rows
               for (i = 0; i < numberDense; i++) {
                    longDouble * a = denseColumn_ + i * numberRows_;
                    for (int j = 0; j < numberRows_; j++) {
                         if (rowsDropped[j])
                              a[j] = 0.0;
                    }
                    for (i = 0; i < numberRows_; i++) {
                         int iRow = permute_[i];
                         workDouble_[i] = a[iRow];
                    }
                    for (i = 0; i < numberRows_; i++) {
                         CoinWorkDouble value = workDouble_[i];
                         CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
                         CoinBigIndex j;
                         for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
                              int iRow = choleskyRow_[j+offset];
                              workDouble_[iRow] -= sparseFactor_[j] * value;
                         }
                    }
                    for (i = 0; i < numberRows_; i++) {
                         int iRow = permute_[i];
                         a[iRow] = workDouble_[i] * diagonal_[i];
                    }
               }
               dense_->resetRowsDropped();
               longDouble * denseBlob = dense_->aMatrix();
               longDouble * denseDiagonal = dense_->diagonal();
               // Update dense matrix
               for (i = 0; i < numberDense; i++) {
                    const longDouble * a = denseColumn_ + i * numberRows_;
                    // do diagonal
                    CoinWorkDouble value = denseDiagonal[i];
                    const longDouble * b = denseColumn_ + i * numberRows_;
                    for (int k = 0; k < numberRows_; k++)
                         value += a[k] * b[k];
                    denseDiagonal[i] = value;
                    for (int j = i + 1; j < numberDense; j++) {
                         CoinWorkDouble value = 0.0;
                         const longDouble * b = denseColumn_ + j * numberRows_;
                         for (int k = 0; k < numberRows_; k++)
                              value += a[k] * b[k];
                         *denseBlob = value;
                         denseBlob++;
                    }
               }
               // dense cholesky (? long double)
               int * dropped = new int [numberDense];
               dense_->factorizePart2(dropped);
               delete [] dropped;
          }
          // try allowing all every time
          //printf("trying ?\n");
          //for (iRow=0;iRow<numberRows_;iRow++) {
          //rowsDropped[iRow]=0;
          //rowsDropped_[iRow]=0;
          //}
          bool cleanCholesky;
          //if (model_->numberIterations()<20||(model_->numberIterations()&1)==0)
          if (model_->numberIterations() < 2000)
               cleanCholesky = true;
          else
               cleanCholesky = false;
          if (cleanCholesky) {
               //drop fresh makes some formADAT easier
               if (newDropped || numberRowsDropped_) {
                    newDropped = 0;
                    for (int i = 0; i < numberRows_; i++) {
                         char dropped = static_cast<char>(rowsDropped[i]);
                         rowsDropped_[i] = dropped;
                         rowsDropped_[i] = 0;
                         if (dropped == 2) {
                              //dropped this time
                              rowsDropped[newDropped++] = i;
                              rowsDropped_[i] = 0;
                         }
                    }
                    numberRowsDropped_ = newDropped;
                    newDropped = -(2 + newDropped);
               }
          } else {
               if (newDropped) {
                    newDropped = 0;
                    for (int i = 0; i < numberRows_; i++) {
                         char dropped = static_cast<char>(rowsDropped[i]);
                         rowsDropped_[i] = dropped;
                         if (dropped == 2) {
                              //dropped this time
                              rowsDropped[newDropped++] = i;
                              rowsDropped_[i] = 1;
                         }
                    }
               }
               numberRowsDropped_ += newDropped;
               if (numberRowsDropped_ && 0) {
                    std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
                              numberRowsDropped_ << " dropped)";
                    if (newDropped) {
                         std::cout << " ( " << newDropped << " dropped this time)";
                    }
                    std::cout << std::endl;
               }
          }
     } else {
          //KKT
          CoinPackedMatrix * quadratic = NULL;
          ClpQuadraticObjective * quadraticObj =
               (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
          if (quadraticObj)
               quadratic = quadraticObj->quadraticObjective();
          // matrix
          int numberRowsModel = model_->numberRows();
          int numberColumns = model_->numberColumns();
          int numberTotal = numberColumns + numberRowsModel;
          // temp
          bool permuted = false;
          for (iRow = 0; iRow < numberRows_; iRow++) {
               if (permute_[iRow] != iRow) {
                    permuted = true;
                    break;
               }
          }
          // but fake it
          for (iRow = 0; iRow < numberRows_; iRow++) {
               //permute_[iRow]=iRow; // force no permute
               //permuteInverse_[permute_[iRow]]=iRow;
          }
          if (permuted) {
               CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
               delta2 *= delta2;
               // Need to permute - ugly
               if (!quadratic) {
                    for (iRow = 0; iRow < numberRows_; iRow++) {
                         longDouble * put = sparseFactor_ + choleskyStart_[iRow];
                         int * which = choleskyRow_ + indexStart_[iRow];
                         int iOriginalRow = permute_[iRow];
                         if (iOriginalRow < numberColumns) {
                              iColumn = iOriginalRow;
                              CoinWorkDouble value = diagonal[iColumn];
                              if (CoinAbs(value) > 1.0e-100) {
                                   value = 1.0 / value;
                                   largest = CoinMax(largest, CoinAbs(value));
                                   diagonal_[iRow] = -value;
                                   CoinBigIndex start = columnStart[iColumn];
                                   CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                                   for (CoinBigIndex j = start; j < end; j++) {
                                        int kRow = row[j] + numberTotal;
                                        kRow = permuteInverse_[kRow];
                                        if (kRow > iRow) {
                                             work[kRow] = element[j];
                                             largest = CoinMax(largest, CoinAbs(element[j]));
                                        }
                                   }
                              } else {
                                   diagonal_[iRow] = -value;
                              }
                         } else if (iOriginalRow < numberTotal) {
                              CoinWorkDouble value = diagonal[iOriginalRow];
                              if (CoinAbs(value) > 1.0e-100) {
                                   value = 1.0 / value;
                                   largest = CoinMax(largest, CoinAbs(value));
                              } else {
                                   value = 1.0e100;
                              }
                              diagonal_[iRow] = -value;
                              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
                              if (kRow > iRow)
                                   work[kRow] = -1.0;
                         } else {
                              diagonal_[iRow] = delta2;
                              int kRow = iOriginalRow - numberTotal;
                              CoinBigIndex start = rowStart[kRow];
                              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = column[j];
                                   int jNewRow = permuteInverse_[jRow];
                                   if (jNewRow > iRow) {
                                        work[jNewRow] = elementByRow[j];
                                        largest = CoinMax(largest, CoinAbs(elementByRow[j]));
                                   }
                              }
                              // slack - should it be permute
                              kRow = permuteInverse_[kRow+numberColumns];
                              if (kRow > iRow)
                                   work[kRow] = -1.0;
                         }
                         CoinBigIndex j;
                         int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
                         for (j = 0; j < number; j++) {
                              int jRow = which[j];
                              put[j] = work[jRow];
                              work[jRow] = 0.0;
                         }
                    }
               } else {
                    // quadratic
                    const int * columnQuadratic = quadratic->getIndices();
                    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
                    const int * columnQuadraticLength = quadratic->getVectorLengths();
                    const double * quadraticElement = quadratic->getElements();
                    for (iRow = 0; iRow < numberRows_; iRow++) {
                         longDouble * put = sparseFactor_ + choleskyStart_[iRow];
                         int * which = choleskyRow_ + indexStart_[iRow];
                         int iOriginalRow = permute_[iRow];
                         if (iOriginalRow < numberColumns) {
                              CoinBigIndex j;
                              iColumn = iOriginalRow;
                              CoinWorkDouble value = diagonal[iColumn];
                              if (CoinAbs(value) > 1.0e-100) {
                                   value = 1.0 / value;
                                   for (j = columnQuadraticStart[iColumn];
                                             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
                                        int jColumn = columnQuadratic[j];
                                        int jNewColumn = permuteInverse_[jColumn];
                                        if (jNewColumn > iRow) {
                                             work[jNewColumn] = -quadraticElement[j];
                                        } else if (iColumn == jColumn) {
                                             value += quadraticElement[j];
                                        }
                                   }
                                   largest = CoinMax(largest, CoinAbs(value));
                                   diagonal_[iRow] = -value;
                                   CoinBigIndex start = columnStart[iColumn];
                                   CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                                   for (j = start; j < end; j++) {
                                        int kRow = row[j] + numberTotal;
                                        kRow = permuteInverse_[kRow];
                                        if (kRow > iRow) {
                                             work[kRow] = element[j];
                                             largest = CoinMax(largest, CoinAbs(element[j]));
                                        }
                                   }
                              } else {
                                   diagonal_[iRow] = -value;
                              }
                         } else if (iOriginalRow < numberTotal) {
                              CoinWorkDouble value = diagonal[iOriginalRow];
                              if (CoinAbs(value) > 1.0e-100) {
                                   value = 1.0 / value;
                                   largest = CoinMax(largest, CoinAbs(value));
                              } else {
                                   value = 1.0e100;
                              }
                              diagonal_[iRow] = -value;
                              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
                              if (kRow > iRow)
                                   work[kRow] = -1.0;
                         } else {
                              diagonal_[iRow] = delta2;
                              int kRow = iOriginalRow - numberTotal;
                              CoinBigIndex start = rowStart[kRow];
                              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   int jRow = column[j];
                                   int jNewRow = permuteInverse_[jRow];
                                   if (jNewRow > iRow) {
                                        work[jNewRow] = elementByRow[j];
                                        largest = CoinMax(largest, CoinAbs(elementByRow[j]));
                                   }
                              }
                              // slack - should it be permute
                              kRow = permuteInverse_[kRow+numberColumns];
                              if (kRow > iRow)
                                   work[kRow] = -1.0;
                         }
                         CoinBigIndex j;
                         int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
                         for (j = 0; j < number; j++) {
                              int jRow = which[j];
                              put[j] = work[jRow];
                              work[jRow] = 0.0;
                         }
                         for (j = 0; j < numberRows_; j++)
                              assert (!work[j]);
                    }
               }
          } else {
               if (!quadratic) {
                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                         longDouble * put = sparseFactor_ + choleskyStart_[iColumn];
                         int * which = choleskyRow_ + indexStart_[iColumn];
                         CoinWorkDouble value = diagonal[iColumn];
                         if (CoinAbs(value) > 1.0e-100) {
                              value = 1.0 / value;
                              largest = CoinMax(largest, CoinAbs(value));
                              diagonal_[iColumn] = -value;
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (CoinBigIndex j = start; j < end; j++) {
                                   //choleskyRow_[numberElements]=row[j]+numberTotal;
                                   //sparseFactor_[numberElements++]=element[j];
                                   work[row[j] + numberTotal] = element[j];
                                   largest = CoinMax(largest, CoinAbs(element[j]));
                              }
                         } else {
                              diagonal_[iColumn] = -value;
                         }
                         CoinBigIndex j;
                         int number = choleskyStart_[iColumn+1] - choleskyStart_[iColumn];
                         for (j = 0; j < number; j++) {
                              int jRow = which[j];
                              put[j] = work[jRow];
                              work[jRow] = 0.0;
                         }
                    }
               } else {
                    // Quadratic
                    const int * columnQuadratic = quadratic->getIndices();
                    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
                    const int * columnQuadraticLength = quadratic->getVectorLengths();
                    const double * quadraticElement = quadratic->getElements();
                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                         longDouble * put = sparseFactor_ + choleskyStart_[iColumn];
                         int * which = choleskyRow_ + indexStart_[iColumn];
                         int number = choleskyStart_[iColumn+1] - choleskyStart_[iColumn];
                         CoinWorkDouble value = diagonal[iColumn];
                         CoinBigIndex j;
                         if (CoinAbs(value) > 1.0e-100) {
                              value = 1.0 / value;
                              for (j = columnQuadraticStart[iColumn];
                                        j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
                                   int jColumn = columnQuadratic[j];
                                   if (jColumn > iColumn) {
                                        work[jColumn] = -quadraticElement[j];
                                   } else if (iColumn == jColumn) {
                                        value += quadraticElement[j];
                                   }
                              }
                              largest = CoinMax(largest, CoinAbs(value));
                              diagonal_[iColumn] = -value;
                              CoinBigIndex start = columnStart[iColumn];
                              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
                              for (j = start; j < end; j++) {
                                   work[row[j] + numberTotal] = element[j];
                                   largest = CoinMax(largest, CoinAbs(element[j]));
                              }
                              for (j = 0; j < number; j++) {
                                   int jRow = which[j];
                                   put[j] = work[jRow];
                                   work[jRow] = 0.0;
                              }
                         } else {
                              value = 1.0e100;
                              diagonal_[iColumn] = -value;
                              for (j = 0; j < number; j++) {
                                   int jRow = which[j];
                                   put[j] = work[jRow];
                              }
                         }
                    }
               }
               // slacks
               for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) {
                    longDouble * put = sparseFactor_ + choleskyStart_[iColumn];
                    int * which = choleskyRow_ + indexStart_[iColumn];
                    CoinWorkDouble value = diagonal[iColumn];
                    if (CoinAbs(value) > 1.0e-100) {
                         value = 1.0 / value;
                         largest = CoinMax(largest, CoinAbs(value));
                    } else {
                         value = 1.0e100;
                    }
                    diagonal_[iColumn] = -value;
                    work[iColumn-numberColumns+numberTotal] = -1.0;
                    CoinBigIndex j;
                    int number = choleskyStart_[iColumn+1] - choleskyStart_[iColumn];
                    for (j = 0; j < number; j++) {
                         int jRow = which[j];
                         put[j] = work[jRow];
                         work[jRow] = 0.0;
                    }
               }
               // Finish diagonal
               CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
               delta2 *= delta2;
               for (iRow = 0; iRow < numberRowsModel; iRow++) {
                    longDouble * put = sparseFactor_ + choleskyStart_[iRow+numberTotal];
                    diagonal_[iRow+numberTotal] = delta2;
                    CoinBigIndex j;
                    int number = choleskyStart_[iRow+numberTotal+1] - choleskyStart_[iRow+numberTotal];
                    for (j = 0; j < number; j++) {
                         put[j] = 0.0;
                    }
               }
          }
          //check sizes
          largest *= 1.0e-20;
          largest = CoinMin(largest, CHOL_SMALL_VALUE);
          doubleParameters_[10] = CoinMax(1.0e-20, largest);
          integerParameters_[20] = 0;
          doubleParameters_[3] = 0.0;
          doubleParameters_[4] = COIN_DBL_MAX;
          // Set up LDL cutoff
          integerParameters_[34] = numberTotal;
          // KKT
          int * rowsDropped2 = new int[numberRows_];
          CoinZeroN(rowsDropped2, numberRows_);
          factorizePart2(rowsDropped2);
          newDropped = integerParameters_[20];
          largest = doubleParameters_[3];
          smallest = doubleParameters_[4];
          if (model_->messageHandler()->logLevel() > 1)
               std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
          choleskyCondition_ = largest / smallest;
          // Should save adjustments in ..R_
          int n1 = 0, n2 = 0;
          CoinWorkDouble * primalR = model_->primalR();
          CoinWorkDouble * dualR = model_->dualR();
          for (iRow = 0; iRow < numberTotal; iRow++) {
               if (rowsDropped2[iRow]) {
                    n1++;
                    //printf("row region1 %d dropped\n",iRow);
                    //rowsDropped_[iRow]=1;
                    rowsDropped_[iRow] = 0;
                    primalR[iRow] = doubleParameters_[20];
               } else {
                    rowsDropped_[iRow] = 0;
                    primalR[iRow] = 0.0;
               }
          }
          for (; iRow < numberRows_; iRow++) {
               if (rowsDropped2[iRow]) {
                    n2++;
                    //printf("row region2 %d dropped\n",iRow);
                    //rowsDropped_[iRow]=1;
                    rowsDropped_[iRow] = 0;
                    dualR[iRow-numberTotal] = doubleParameters_[34];
               } else {
                    rowsDropped_[iRow] = 0;
                    dualR[iRow-numberTotal] = 0.0;
               }
          }
          delete [] rowsDropped2;
     }
     status_ = 0;
     return newDropped;
}
/* Factorize - filling in rowsDropped and returning number dropped
   in integerParam.
*/
void
ClpCholeskyBase::factorizePart2(int * rowsDropped)
{
     CoinWorkDouble largest = doubleParameters_[3];
     CoinWorkDouble smallest = doubleParameters_[4];
     // probably done before
     largest = 0.0;
     smallest = COIN_DBL_MAX;
     double dropValue = doubleParameters_[10];
     int firstPositive = integerParameters_[34];
     longDouble * d = ClpCopyOfArray(diagonal_, numberRows_);
     int iRow;
     // minimum size before clique done
     //#define MINCLIQUE INT_MAX
#define MINCLIQUE 3
     longDouble * work = workDouble_;
     CoinBigIndex * first = workInteger_;

     for (iRow = 0; iRow < numberRows_; iRow++) {
          link_[iRow] = -1;
          work[iRow] = 0.0;
          first[iRow] = choleskyStart_[iRow];
     }

     int lastClique = -1;
     bool inClique = false;
     bool newClique = false;
     bool endClique = false;
     int lastRow = 0;
     int nextRow2 = -1;

     for (iRow = 0; iRow < firstDense_ + 1; iRow++) {
          if (iRow < firstDense_) {
               endClique = false;
               if (clique_[iRow] > 0) {
                    // this is in a clique
                    inClique = true;
                    if (clique_[iRow] > lastClique) {
                         // new Clique
                         newClique = true;
                         // If we have clique going then signal to do old one
                         endClique = (lastClique > 0);
                    } else {
                         // Still in clique
                         newClique = false;
                    }
               } else {
                    // not in clique
                    inClique = false;
                    newClique = false;
                    // If we have clique going then signal to do old one
                    endClique = (lastClique > 0);
               }
               lastClique = clique_[iRow];
          } else if (inClique) {
               // Finish off
               endClique = true;
          } else {
               break;
          }
          if (endClique) {
               // We have just finished updating a clique - do block pivot and clean up
               int jRow;
               for ( jRow = lastRow; jRow < iRow; jRow++) {
                    int jCount = jRow - lastRow;
                    CoinWorkDouble diagonalValue = diagonal_[jRow];
                    CoinBigIndex start = choleskyStart_[jRow];
                    CoinBigIndex end = choleskyStart_[jRow+1];
                    for (int kRow = lastRow; kRow < jRow; kRow++) {
                         jCount--;
                         CoinBigIndex get = choleskyStart_[kRow] + jCount;
                         CoinWorkDouble a_jk = sparseFactor_[get];
                         CoinWorkDouble value1 = d[kRow] * a_jk;
                         diagonalValue -= a_jk * value1;
                         for (CoinBigIndex j = start; j < end; j++)
                              sparseFactor_[j] -= value1 * sparseFactor_[++get];
                    }
                    // check
                    int originalRow = permute_[jRow];
                    if (originalRow < firstPositive) {
                         // must be negative
                         if (diagonalValue <= -dropValue) {
                              smallest = CoinMin(smallest, -diagonalValue);
                              largest = CoinMax(largest, -diagonalValue);
                              d[jRow] = diagonalValue;
                              diagonalValue = 1.0 / diagonalValue;
                         } else {
                              rowsDropped[originalRow] = 2;
                              d[jRow] = -1.0e100;
                              diagonalValue = 0.0;
                              integerParameters_[20]++;
                         }
                    } else {
                         // must be positive
                         if (diagonalValue >= dropValue) {
                              smallest = CoinMin(smallest, diagonalValue);
                              largest = CoinMax(largest, diagonalValue);
                              d[jRow] = diagonalValue;
                              diagonalValue = 1.0 / diagonalValue;
                         } else {
                              rowsDropped[originalRow] = 2;
                              d[jRow] = 1.0e100;
                              diagonalValue = 0.0;
                              integerParameters_[20]++;
                         }
                    }
                    diagonal_[jRow] = diagonalValue;
                    for (CoinBigIndex j = start; j < end; j++) {
                         sparseFactor_[j] *= diagonalValue;
                    }
               }
               if (nextRow2 >= 0) {
                    for ( jRow = lastRow; jRow < iRow - 1; jRow++) {
                         link_[jRow] = jRow + 1;
                    }
                    link_[iRow-1] = link_[nextRow2];
                    link_[nextRow2] = lastRow;
               }
          }
          if (iRow == firstDense_)
               break; // we were just cleaning up
          if (newClique) {
               // initialize new clique
               lastRow = iRow;
          }
          // for each column L[*,kRow] that affects L[*,iRow]
          CoinWorkDouble diagonalValue = diagonal_[iRow];
          int nextRow = link_[iRow];
          int kRow = 0;
          while (1) {
               kRow = nextRow;
               if (kRow < 0)
                    break; // out of loop
               nextRow = link_[kRow];
               // Modify by outer product of L[*,irow] by L[*,krow] from first
               CoinBigIndex k = first[kRow];
               CoinBigIndex end = choleskyStart_[kRow+1];
               assert(k < end);
               CoinWorkDouble a_ik = sparseFactor_[k++];
               CoinWorkDouble value1 = d[kRow] * a_ik;
               // update first
               first[kRow] = k;
               diagonalValue -= value1 * a_ik;
               CoinBigIndex offset = indexStart_[kRow] - choleskyStart_[kRow];
               if (k < end) {
                    int jRow = choleskyRow_[k+offset];
                    if (clique_[kRow] < MINCLIQUE) {
                         link_[kRow] = link_[jRow];
                         link_[jRow] = kRow;
                         for (; k < end; k++) {
                              int jRow = choleskyRow_[k+offset];
                              work[jRow] += sparseFactor_[k] * value1;
                         }
                    } else {
                         // Clique
                         CoinBigIndex currentIndex = k + offset;
                         int linkSave = link_[jRow];
                         link_[jRow] = kRow;
                         work[kRow] = value1; // ? or a_jk
                         int last = kRow + clique_[kRow];
                         for (int kkRow = kRow + 1; kkRow < last; kkRow++) {
                              CoinBigIndex j = first[kkRow];
                              //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
                              CoinWorkDouble a = sparseFactor_[j];
                              CoinWorkDouble dValue = d[kkRow] * a;
                              diagonalValue -= a * dValue;
                              work[kkRow] = dValue;
                              first[kkRow]++;
                              link_[kkRow-1] = kkRow;
                         }
                         nextRow = link_[last-1];
                         link_[last-1] = linkSave;
                         int length = end - k;
                         for (int i = 0; i < length; i++) {
                              int lRow = choleskyRow_[currentIndex++];
                              CoinWorkDouble t0 = work[lRow];
                              for (int kkRow = kRow; kkRow < last; kkRow++) {
                                   CoinBigIndex j = first[kkRow] + i;
                                   t0 += work[kkRow] * sparseFactor_[j];
                              }
                              work[lRow] = t0;
                         }
                    }
               }
          }
          // Now apply
          if (inClique) {
               // in clique
               diagonal_[iRow] = diagonalValue;
               CoinBigIndex start = choleskyStart_[iRow];
               CoinBigIndex end = choleskyStart_[iRow+1];
               CoinBigIndex currentIndex = indexStart_[iRow];
               nextRow2 = -1;
               CoinBigIndex get = start + clique_[iRow] - 1;
               if (get < end) {
                    nextRow2 = choleskyRow_[currentIndex+get-start];
                    first[iRow] = get;
               }
               for (CoinBigIndex j = start; j < end; j++) {
                    int kRow = choleskyRow_[currentIndex++];
                    sparseFactor_[j] -= work[kRow]; // times?
                    work[kRow] = 0.0;
               }
          } else {
               // not in clique
               int originalRow = permute_[iRow];
               if (originalRow < firstPositive) {
                    // must be negative
                    if (diagonalValue <= -dropValue) {
                         smallest = CoinMin(smallest, -diagonalValue);
                         largest = CoinMax(largest, -diagonalValue);
                         d[iRow] = diagonalValue;
                         diagonalValue = 1.0 / diagonalValue;
                    } else {
                         rowsDropped[originalRow] = 2;
                         d[iRow] = -1.0e100;
                         diagonalValue = 0.0;
                         integerParameters_[20]++;
                    }
               } else {
                    // must be positive
                    if (diagonalValue >= dropValue) {
                         smallest = CoinMin(smallest, diagonalValue);
                         largest = CoinMax(largest, diagonalValue);
                         d[iRow] = diagonalValue;
                         diagonalValue = 1.0 / diagonalValue;
                    } else {
                         rowsDropped[originalRow] = 2;
                         d[iRow] = 1.0e100;
                         diagonalValue = 0.0;
                         integerParameters_[20]++;
                    }
               }
               diagonal_[iRow] = diagonalValue;
               CoinBigIndex offset = indexStart_[iRow] - choleskyStart_[iRow];
               CoinBigIndex start = choleskyStart_[iRow];
               CoinBigIndex end = choleskyStart_[iRow+1];
               assert (first[iRow] == start);
               if (start < end) {
                    int nextRow = choleskyRow_[start+offset];
                    link_[iRow] = link_[nextRow];
                    link_[nextRow] = iRow;
                    for (CoinBigIndex j = start; j < end; j++) {
                         int jRow = choleskyRow_[j+offset];
                         CoinWorkDouble value = sparseFactor_[j] - work[jRow];
                         work[jRow] = 0.0;
                         sparseFactor_[j] = diagonalValue * value;
                    }
               }
          }
     }
     if (firstDense_ < numberRows_) {
          // do dense
          // update dense part
          updateDense(d,/*work,*/first);
          ClpCholeskyDense dense;
          // just borrow space
          int nDense = numberRows_ - firstDense_;
          if (doKKT_) {
               for (iRow = firstDense_; iRow < numberRows_; iRow++) {
                    int originalRow = permute_[iRow];
                    if (originalRow >= firstPositive) {
                         firstPositive = iRow - firstDense_;
                         break;
                    }
               }
          }
          dense.reserveSpace(this, nDense);
          int * dropped = new int[nDense];
          memset(dropped, 0, nDense * sizeof(int));
          dense.setDoubleParameter(3, largest);
          dense.setDoubleParameter(4, smallest);
          dense.setDoubleParameter(10, dropValue);
          dense.setIntegerParameter(20, 0);
          dense.setIntegerParameter(34, firstPositive);
          dense.setModel(model_);
          dense.factorizePart2(dropped);
          largest = dense.getDoubleParameter(3);
          smallest = dense.getDoubleParameter(4);
          integerParameters_[20] += dense.getIntegerParameter(20);
          for (iRow = firstDense_; iRow < numberRows_; iRow++) {
               int originalRow = permute_[iRow];
               rowsDropped[originalRow] = dropped[iRow-firstDense_];
          }
          delete [] dropped;
     }
     delete [] d;
     doubleParameters_[3] = largest;
     doubleParameters_[4] = smallest;
     return;
}
// Updates dense part (broken out for profiling)
void ClpCholeskyBase::updateDense(longDouble * d, /*longDouble * work,*/ int * first)
{
     for (int iRow = 0; iRow < firstDense_; iRow++) {
          CoinBigIndex start = first[iRow];
          CoinBigIndex end = choleskyStart_[iRow+1];
          if (start < end) {
               CoinBigIndex offset = indexStart_[iRow] - choleskyStart_[iRow];
               if (clique_[iRow] < 2) {
                    CoinWorkDouble dValue = d[iRow];
                    for (CoinBigIndex k = start; k < end; k++) {
                         int kRow = choleskyRow_[k+offset];
                         assert(kRow >= firstDense_);
                         CoinWorkDouble a_ik = sparseFactor_[k];
                         CoinWorkDouble value1 = dValue * a_ik;
                         diagonal_[kRow] -= value1 * a_ik;
                         CoinBigIndex base = choleskyStart_[kRow] - kRow - 1;
                         for (CoinBigIndex j = k + 1; j < end; j++) {
                              int jRow = choleskyRow_[j+offset];
                              CoinWorkDouble a_jk = sparseFactor_[j];
                              sparseFactor_[base+jRow] -= a_jk * value1;
                         }
                    }
               } else if (clique_[iRow] < 3) {
                    // do as pair
                    CoinWorkDouble dValue0 = d[iRow];
                    CoinWorkDouble dValue1 = d[iRow+1];
                    int offset1 = first[iRow+1] - start;
                    // skip row
                    iRow++;
                    for (CoinBigIndex k = start; k < end; k++) {
                         int kRow = choleskyRow_[k+offset];
                         assert(kRow >= firstDense_);
                         CoinWorkDouble a_ik0 = sparseFactor_[k];
                         CoinWorkDouble value0 = dValue0 * a_ik0;
                         CoinWorkDouble a_ik1 = sparseFactor_[k+offset1];
                         CoinWorkDouble value1 = dValue1 * a_ik1;
                         diagonal_[kRow] -= value0 * a_ik0 + value1 * a_ik1;
                         CoinBigIndex base = choleskyStart_[kRow] - kRow - 1;
                         for (CoinBigIndex j = k + 1; j < end; j++) {
                              int jRow = choleskyRow_[j+offset];
                              CoinWorkDouble a_jk0 = sparseFactor_[j];
                              CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
                              sparseFactor_[base+jRow] -= a_jk0 * value0 + a_jk1 * value1;
                         }
                    }
#define MANY_REGISTERS
#ifdef MANY_REGISTERS
               } else if (clique_[iRow] == 3) {
#else
               } else {
#endif
                    // do as clique
                    // maybe later get fancy on big cliques and do transpose copy
                    // seems only worth going to 3 on Intel
                    CoinWorkDouble dValue0 = d[iRow];
                    CoinWorkDouble dValue1 = d[iRow+1];
                    CoinWorkDouble dValue2 = d[iRow+2];
                    // get offsets and skip rows
                    int offset1 = first[++iRow] - start;
                    int offset2 = first[++iRow] - start;
                    for (CoinBigIndex k = start; k < end; k++) {
                         int kRow = choleskyRow_[k+offset];
                         assert(kRow >= firstDense_);
                         CoinWorkDouble diagonalValue = diagonal_[kRow];
                         CoinWorkDouble a_ik0 = sparseFactor_[k];
                         CoinWorkDouble value0 = dValue0 * a_ik0;
                         CoinWorkDouble a_ik1 = sparseFactor_[k+offset1];
                         CoinWorkDouble value1 = dValue1 * a_ik1;
                         CoinWorkDouble a_ik2 = sparseFactor_[k+offset2];
                         CoinWorkDouble value2 = dValue2 * a_ik2;
                         CoinBigIndex base = choleskyStart_[kRow] - kRow - 1;
                         diagonal_[kRow] = diagonalValue - value0 * a_ik0 - value1 * a_ik1 - value2 * a_ik2;
                         for (CoinBigIndex j = k + 1; j < end; j++) {
                              int jRow = choleskyRow_[j+offset];
                              CoinWorkDouble a_jk0 = sparseFactor_[j];
                              CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
                              CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
                              sparseFactor_[base+jRow] -= a_jk0 * value0 + a_jk1 * value1 + a_jk2 * value2;
                         }
                    }
#ifdef MANY_REGISTERS
               }
               else {
                    // do as clique
                    // maybe later get fancy on big cliques and do transpose copy
                    // maybe only worth going to 3 on Intel (but may have hidden registers)
                    CoinWorkDouble dValue0 = d[iRow];
                    CoinWorkDouble dValue1 = d[iRow+1];
                    CoinWorkDouble dValue2 = d[iRow+2];
                    CoinWorkDouble dValue3 = d[iRow+3];
                    // get offsets and skip rows
                    int offset1 = first[++iRow] - start;
                    int offset2 = first[++iRow] - start;
                    int offset3 = first[++iRow] - start;
                    for (CoinBigIndex k = start; k < end; k++) {
                         int kRow = choleskyRow_[k+offset];
                         assert(kRow >= firstDense_);
                         CoinWorkDouble diagonalValue = diagonal_[kRow];
                         CoinWorkDouble a_ik0 = sparseFactor_[k];
                         CoinWorkDouble value0 = dValue0 * a_ik0;
                         CoinWorkDouble a_ik1 = sparseFactor_[k+offset1];
                         CoinWorkDouble value1 = dValue1 * a_ik1;
                         CoinWorkDouble a_ik2 = sparseFactor_[k+offset2];
                         CoinWorkDouble value2 = dValue2 * a_ik2;
                         CoinWorkDouble a_ik3 = sparseFactor_[k+offset3];
                         CoinWorkDouble value3 = dValue3 * a_ik3;
                         CoinBigIndex base = choleskyStart_[kRow] - kRow - 1;
                         diagonal_[kRow] = diagonalValue - (value0 * a_ik0 + value1 * a_ik1 + value2 * a_ik2 + value3 * a_ik3);
                         for (CoinBigIndex j = k + 1; j < end; j++) {
                              int jRow = choleskyRow_[j+offset];
                              CoinWorkDouble a_jk0 = sparseFactor_[j];
                              CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
                              CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
                              CoinWorkDouble a_jk3 = sparseFactor_[j+offset3];
                              sparseFactor_[base+jRow] -= a_jk0 * value0 + a_jk1 * value1 + a_jk2 * value2 + a_jk3 * value3;
                         }
                    }
#endif
               }
          }
     }
}
/* Uses factorization to solve. */
void
ClpCholeskyBase::solve (CoinWorkDouble * region)
{
     if (!whichDense_) {
          solve(region, 3);
     } else {
          // dense columns
          int i;
          solve(region, 1);
          // do change;
          int numberDense = dense_->numberRows();
          CoinWorkDouble * change = new CoinWorkDouble[numberDense];
          for (i = 0; i < numberDense; i++) {
               const longDouble * a = denseColumn_ + i * numberRows_;
               CoinWorkDouble value = 0.0;
               for (int iRow = 0; iRow < numberRows_; iRow++)
                    value += a[iRow] * region[iRow];
               change[i] = value;
          }
          // solve
          dense_->solve(change);
          for (i = 0; i < numberDense; i++) {
               const longDouble * a = denseColumn_ + i * numberRows_;
               CoinWorkDouble value = change[i];
               for (int iRow = 0; iRow < numberRows_; iRow++)
                    region[iRow] -= value * a[iRow];
          }
          delete [] change;
          // and finish off
          solve(region, 2);
     }
}
/* solve - 1 just first half, 2 just second half - 3 both.
   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
*/
void
ClpCholeskyBase::solve(CoinWorkDouble * region, int type)
{
#ifdef CLP_DEBUG
     double * regionX = NULL;
     if (sizeof(CoinWorkDouble) != sizeof(double) && type == 3) {
          regionX = ClpCopyOfArray(region, numberRows_);
     }
#endif
     CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
     int i;
     CoinBigIndex j;
     for (i = 0; i < numberRows_; i++) {
          int iRow = permute_[i];
          work[i] = region[iRow];
     }
     switch (type) {
     case 1:
          for (i = 0; i < numberRows_; i++) {
               CoinWorkDouble value = work[i];
               CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
               for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
                    int iRow = choleskyRow_[j+offset];
                    work[iRow] -= sparseFactor_[j] * value;
               }
          }
          for (i = 0; i < numberRows_; i++) {
               int iRow = permute_[i];
               region[iRow] = work[i] * diagonal_[i];
          }
          break;
     case 2:
          for (i = numberRows_ - 1; i >= 0; i--) {
               CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
               CoinWorkDouble value = work[i] * diagonal_[i];
               for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
                    int iRow = choleskyRow_[j+offset];
                    value -= sparseFactor_[j] * work[iRow];
               }
               work[i] = value;
               int iRow = permute_[i];
               region[iRow] = value;
          }
          break;
     case 3:
          for (i = 0; i < firstDense_; i++) {
               CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
               CoinWorkDouble value = work[i];
               for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
                    int iRow = choleskyRow_[j+offset];
                    work[iRow] -= sparseFactor_[j] * value;
               }
          }
          if (firstDense_ < numberRows_) {
               // do dense
               ClpCholeskyDense dense;
               // just borrow space
               int nDense = numberRows_ - firstDense_;
               dense.reserveSpace(this, nDense);
               dense.solve(work + firstDense_);
               for (i = numberRows_ - 1; i >= firstDense_; i--) {
                    CoinWorkDouble value = work[i];
                    int iRow = permute_[i];
                    region[iRow] = value;
               }
          }
          for (i = firstDense_ - 1; i >= 0; i--) {
               CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
               CoinWorkDouble value = work[i] * diagonal_[i];
               for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
                    int iRow = choleskyRow_[j+offset];
                    value -= sparseFactor_[j] * work[iRow];
               }
               work[i] = value;
               int iRow = permute_[i];
               region[iRow] = value;
          }
          break;
     }
#ifdef CLP_DEBUG
     if (regionX) {
          longDouble * work = workDouble_;
          int i;
          CoinBigIndex j;
          double largestO = 0.0;
          for (i = 0; i < numberRows_; i++) {
               largestO = CoinMax(largestO, CoinAbs(regionX[i]));
          }
          for (i = 0; i < numberRows_; i++) {
               int iRow = permute_[i];
               work[i] = regionX[iRow];
          }
          for (i = 0; i < firstDense_; i++) {
               CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
               CoinWorkDouble value = work[i];
               for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
                    int iRow = choleskyRow_[j+offset];
                    work[iRow] -= sparseFactor_[j] * value;
               }
          }
          if (firstDense_ < numberRows_) {
               // do dense
               ClpCholeskyDense dense;
               // just borrow space
               int nDense = numberRows_ - firstDense_;
               dense.reserveSpace(this, nDense);
               dense.solve(work + firstDense_);
               for (i = numberRows_ - 1; i >= firstDense_; i--) {
                    CoinWorkDouble value = work[i];
                    int iRow = permute_[i];
                    regionX[iRow] = value;
               }
          }
          for (i = firstDense_ - 1; i >= 0; i--) {
               CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
               CoinWorkDouble value = work[i] * diagonal_[i];
               for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
                    int iRow = choleskyRow_[j+offset];
                    value -= sparseFactor_[j] * work[iRow];
               }
               work[i] = value;
               int iRow = permute_[i];
               regionX[iRow] = value;
          }
          double largest = 0.0;
          double largestV = 0.0;
          for (i = 0; i < numberRows_; i++) {
               largest = CoinMax(largest, CoinAbs(region[i] - regionX[i]));
               largestV = CoinMax(largestV, CoinAbs(region[i]));
          }
          printf("largest difference %g, largest %g, largest original %g\n",
                 largest, largestV, largestO);
          delete [] regionX;
     }
#endif
}
#if 0 //CLP_LONG_CHOLESKY
/* Uses factorization to solve. */
void
ClpCholeskyBase::solve (CoinWorkDouble * region)
{
     assert (!whichDense_) ;
     CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
     int i;
     CoinBigIndex j;
     for (i = 0; i < numberRows_; i++) {
          int iRow = permute_[i];
          work[i] = region[iRow];
     }
     for (i = 0; i < firstDense_; i++) {
          CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
          CoinWorkDouble value = work[i];
          for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
               int iRow = choleskyRow_[j+offset];
               work[iRow] -= sparseFactor_[j] * value;
          }
     }
     if (firstDense_ < numberRows_) {
          // do dense
          ClpCholeskyDense dense;
          // just borrow space
          int nDense = numberRows_ - firstDense_;
          dense.reserveSpace(this, nDense);
          dense.solve(work + firstDense_);
          for (i = numberRows_ - 1; i >= firstDense_; i--) {
               CoinWorkDouble value = work[i];
               int iRow = permute_[i];
               region[iRow] = value;
          }
     }
     for (i = firstDense_ - 1; i >= 0; i--) {
          CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
          CoinWorkDouble value = work[i] * diagonal_[i];
          for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
               int iRow = choleskyRow_[j+offset];
               value -= sparseFactor_[j] * work[iRow];
          }
          work[i] = value;
          int iRow = permute_[i];
          region[iRow] = value;
     }
}
#endif