Changeset 414


Ignore:
Timestamp:
Aug 29, 2004 6:09:18 AM (15 years ago)
Author:
forrest
Message:

barrier stuff

Location:
trunk
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpCholeskyBase.cpp

    r328 r414  
    1010#include "ClpHelperFunctions.hpp"
    1111#include "CoinHelperFunctions.hpp"
     12#include "CoinSort.hpp"
     13#include "ClpCholeskyDense.hpp"
     14#include "ClpMessage.hpp"
     15#include "ClpQuadraticObjective.hpp"
    1216
    1317//#############################################################################
     
    1822// Default Constructor
    1923//-------------------------------------------------------------------
    20 ClpCholeskyBase::ClpCholeskyBase () :
    21   type_(-1),
    22   pivotTolerance_(1.0e-14),
     24ClpCholeskyBase::ClpCholeskyBase (int denseThreshold) :
     25  type_(0),
     26  doKKT_(false),
    2327  zeroTolerance_(1.0e-17),
    2428  choleskyCondition_(0.0),
     
    2832  status_(0),
    2933  rowsDropped_(NULL),
    30   permuteIn_(NULL),
    31   permuteOut_(NULL),
    32   numberRowsDropped_(0)
     34  permuteInverse_(NULL),
     35  permute_(NULL),
     36  numberRowsDropped_(0),
     37  sparseFactor_(NULL),
     38  choleskyStart_(NULL),
     39  choleskyRow_(NULL),
     40  indexStart_(NULL),
     41  diagonal_(NULL),
     42  workDouble_(NULL),
     43  link_(NULL),
     44  workInteger_(NULL),
     45  clique_(NULL),
     46  sizeFactor_(0),
     47  sizeIndex_(0),
     48  firstDense_(0),
     49  rowCopy_(NULL),
     50  whichDense_(NULL),
     51  denseColumn_(NULL),
     52  dense_(NULL),
     53  denseThreshold_(denseThreshold)
    3354{
    34 
     55  memset(integerParameters_,0,64*sizeof(int));
     56  memset(doubleParameters_,0,64*sizeof(double));
    3557}
    3658
     
    4062ClpCholeskyBase::ClpCholeskyBase (const ClpCholeskyBase & rhs) :
    4163  type_(rhs.type_),
    42   pivotTolerance_(rhs.pivotTolerance_),
     64  doKKT_(rhs.doKKT_),
    4365  zeroTolerance_(rhs.zeroTolerance_),
    4466  choleskyCondition_(rhs.choleskyCondition_),
     
    5072
    5173  rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_,numberRows_);
    52   permuteIn_ = ClpCopyOfArray(rhs.permuteIn_,numberRows_);
    53   permuteOut_ = ClpCopyOfArray(rhs.permuteOut_,numberRows_);
     74  permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_,numberRows_);
     75  permute_ = ClpCopyOfArray(rhs.permute_,numberRows_);
     76  sizeFactor_=rhs.sizeFactor_;
     77  sizeIndex_ = rhs.sizeIndex_;
     78  firstDense_ = rhs.firstDense_;
     79  sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
     80  choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
     81  indexStart_ = ClpCopyOfArray(rhs.indexStart_,numberRows_);
     82  choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
     83  diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
     84  workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
     85  link_ = ClpCopyOfArray(rhs.link_,numberRows_);
     86  workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
     87  clique_ = ClpCopyOfArray(rhs.clique_,numberRows_);
     88  memcpy(integerParameters_,rhs.integerParameters_,64*sizeof(int));
     89  memcpy(doubleParameters_,rhs.doubleParameters_,64*sizeof(double));
     90  rowCopy_ = rhs.rowCopy_->clone();
     91  whichDense_ = NULL;
     92  denseColumn_=NULL;
     93  dense_=NULL;
     94  denseThreshold_ = rhs.denseThreshold_;
    5495}
    5596
     
    60101{
    61102  delete [] rowsDropped_;
    62   delete [] permuteIn_;
    63   delete [] permuteOut_;
     103  delete [] permuteInverse_;
     104  delete [] permute_;
     105  delete [] sparseFactor_;
     106  delete [] choleskyStart_;
     107  delete [] choleskyRow_;
     108  delete [] indexStart_;
     109  delete [] diagonal_;
     110  delete [] workDouble_;
     111  delete [] link_;
     112  delete [] workInteger_;
     113  delete [] clique_;
     114  delete rowCopy_;
     115  delete [] whichDense_;
     116  delete [] denseColumn_;
     117  delete dense_;
    64118}
    65119
     
    72126  if (this != &rhs) {
    73127    type_ = rhs.type_;
    74     pivotTolerance_ = rhs.pivotTolerance_;
     128    doKKT_ = rhs.doKKT_;
    75129    zeroTolerance_ = rhs.zeroTolerance_;
    76130    choleskyCondition_ = rhs.choleskyCondition_;
     
    81135    numberRowsDropped_ = rhs.numberRowsDropped_;
    82136    delete [] rowsDropped_;
    83     delete [] permuteIn_;
    84     delete [] permuteOut_;
     137    delete [] permuteInverse_;
     138    delete [] permute_;
     139    delete [] sparseFactor_;
     140    delete [] choleskyStart_;
     141    delete [] choleskyRow_;
     142    delete [] indexStart_;
     143    delete [] diagonal_;
     144    delete [] workDouble_;
     145    delete [] link_;
     146    delete [] workInteger_;
     147    delete [] clique_;
     148    delete rowCopy_;
     149    delete [] whichDense_;
     150    delete [] denseColumn_;
     151    delete dense_;
    85152    rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_,numberRows_);
    86     permuteIn_ = ClpCopyOfArray(rhs.permuteIn_,numberRows_);
    87     permuteOut_ = ClpCopyOfArray(rhs.permuteOut_,numberRows_);
     153    permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_,numberRows_);
     154    permute_ = ClpCopyOfArray(rhs.permute_,numberRows_);
     155    sizeFactor_=rhs.sizeFactor_;
     156    sizeIndex_ = rhs.sizeIndex_;
     157    firstDense_ = rhs.firstDense_;
     158    sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
     159    choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
     160    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
     161    indexStart_ = ClpCopyOfArray(rhs.indexStart_,numberRows_);
     162    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
     163    diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
     164    workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
     165    link_ = ClpCopyOfArray(rhs.link_,numberRows_);
     166    workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
     167    clique_ = ClpCopyOfArray(rhs.clique_,numberRows_);
     168    delete rowCopy_;
     169    rowCopy_ = rhs.rowCopy_->clone();
     170    whichDense_ = NULL;
     171    denseColumn_=NULL;
     172    dense_=NULL;
     173    denseThreshold_ = rhs.denseThreshold_;
    88174  }
    89175  return *this;
     
    102188                           double diagonalScaleFactor)
    103189{
     190  if (!doKKT_) {
     191    int iColumn;
     192    int numberColumns = model_->numberColumns();
     193    int numberTotal = numberRows_+numberColumns;
     194    double * region1Save = new double[numberTotal];
     195    for (iColumn=0;iColumn<numberTotal;iColumn++) {
     196      region1[iColumn] *= diagonal[iColumn];
     197      region1Save[iColumn]=region1[iColumn];
     198    }
     199    multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0);
     200    model_->clpMatrix()->times(1.0,region1,region2);
     201    double maximumRHS = maximumAbsElement(region2,numberRows_);
     202    double scale=1.0;
     203    double unscale=1.0;
     204    if (maximumRHS>1.0e-30) {
     205      if (maximumRHS<=0.5) {
     206        double factor=2.0;
     207        while (maximumRHS<=0.5) {
     208          maximumRHS*=factor;
     209          scale*=factor;
     210        } /* endwhile */
     211      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
     212        double factor=0.5;
     213        while (maximumRHS>=2.0) {
     214          maximumRHS*=factor;
     215          scale*=factor;
     216        } /* endwhile */
     217      }
     218      unscale=diagonalScaleFactor/scale;
     219    } else {
     220      //effectively zero
     221      scale=0.0;
     222      unscale=0.0;
     223    }
     224    multiplyAdd(NULL,numberRows_,0.0,region2,scale);
     225    solve(region2);
     226    multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
     227    multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns,0.0);
     228    CoinZeroN(region1,numberColumns);
     229    model_->clpMatrix()->transposeTimes(1.0,region2,region1);
     230    for (iColumn=0;iColumn<numberTotal;iColumn++)
     231      region1[iColumn] = region1[iColumn]*diagonal[iColumn]-region1Save[iColumn];
     232    delete [] region1Save;
     233  } else {
     234    // KKT
     235    int numberRowsModel = model_->numberRows();
     236    int numberColumns = model_->numberColumns();
     237    int numberTotal = numberColumns + numberRowsModel;
     238    double * array = new double [numberRows_];
     239    CoinMemcpyN(region1,numberTotal,array);
     240    CoinMemcpyN(region2,numberRowsModel,array+numberTotal);
     241    solve(array);
     242    int iRow;
     243    for (iRow=0;iRow<numberTotal;iRow++) {
     244      if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
     245        printf("row region1 %d dropped %g\n",iRow,array[iRow]);
     246      }
     247    }
     248    for (;iRow<numberRows_;iRow++) {
     249      if (rowsDropped_[iRow]&&fabs(array[iRow])>1.0e-8) {
     250        printf("row region2 %d dropped %g\n",iRow,array[iRow]);
     251      }
     252    }
     253    CoinMemcpyN(array+numberTotal,numberRowsModel,region2);
     254    CoinMemcpyN(array,numberTotal,region1);
     255    delete [] array;
     256  }
     257}
     258//-------------------------------------------------------------------
     259// Clone
     260//-------------------------------------------------------------------
     261ClpCholeskyBase * ClpCholeskyBase::clone() const
     262{
     263  return new ClpCholeskyBase(*this);
     264}
     265// Forms ADAT - returns nonzero if not enough memory
     266int
     267ClpCholeskyBase::preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT)
     268{
     269  delete rowCopy_;
     270  rowCopy_ = model_->clpMatrix()->reverseOrderedCopy();
     271  if (!doKKT) {
     272    numberRows_ = model_->numberRows();
     273    rowsDropped_ = new char [numberRows_];
     274    memset(rowsDropped_,0,numberRows_);
     275    numberRowsDropped_=0;
     276    // Space for starts
     277    choleskyStart_ = new CoinBigIndex[numberRows_+1];
     278    const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     279    const int * columnLength = model_->clpMatrix()->getVectorLengths();
     280    const int * row = model_->clpMatrix()->getIndices();
     281    const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
     282    const int * rowLength = rowCopy_->getVectorLengths();
     283    const int * column = rowCopy_->getIndices();
     284    // We need two arrays for counts
     285    int * which = new int [numberRows_];
     286    int * used = new int[numberRows_+1];
     287    CoinZeroN(used,numberRows_);
     288    int iRow;
     289    sizeFactor_=0;
     290    int numberColumns = model_->numberColumns();
     291    int numberDense=0;
     292    //denseThreshold_=3;
     293    if (denseThreshold_>0) {
     294      delete [] whichDense_;
     295      delete [] denseColumn_;
     296      delete dense_;
     297      whichDense_ = new char[numberColumns];
     298      int iColumn;
     299      used[numberRows_]=0;
     300      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     301        int length = columnLength[iColumn];
     302        used[length] += 1;
     303      }
     304      int nLong=0;
     305      int stop = CoinMax(denseThreshold_/2,100);
     306      for (iRow=numberRows_;iRow>=stop;iRow--) {
     307        if (used[iRow])
     308          printf("%d columns are of length %d\n",used[iRow],iRow);
     309        nLong += used[iRow];
     310        if (nLong>50||nLong>(numberColumns>>2))
     311          break;
     312      }
     313      CoinZeroN(used,numberRows_);
     314      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     315        if (columnLength[iColumn]<denseThreshold_) {
     316          whichDense_[iColumn]=0;
     317        } else {
     318          whichDense_[iColumn]=1;
     319          numberDense++;
     320        }
     321      }
     322      if (!numberDense||numberDense>100) {
     323        // free
     324        delete [] whichDense_;
     325        whichDense_=NULL;
     326        denseColumn_=NULL;
     327        dense_=NULL;
     328      } else {
     329        // space for dense columns
     330        denseColumn_ = new double [numberDense*numberRows_];
     331        // dense cholesky
     332        dense_ = new ClpCholeskyDense();
     333        dense_->reserveSpace(NULL,numberDense);
     334        printf("Taking %d columns as dense\n",numberDense);
     335      }
     336    }
     337    int offset = includeDiagonal ? 0 : 1;
     338    if (lowerTriangular)
     339      offset=-offset;
     340    for (iRow=0;iRow<numberRows_;iRow++) {
     341      int number=0;
     342      // make sure diagonal exists if includeDiagonal
     343      if (!offset) {
     344        which[0]=iRow;
     345        used[iRow]=1;
     346        number=1;
     347      }
     348      CoinBigIndex startRow=rowStart[iRow];
     349      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
     350      if (lowerTriangular) {
     351        for (CoinBigIndex k=startRow;k<endRow;k++) {
     352          int iColumn=column[k];
     353          if (!whichDense_||!whichDense_[iColumn]) {
     354            CoinBigIndex start=columnStart[iColumn];
     355            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     356            for (CoinBigIndex j=start;j<end;j++) {
     357              int jRow=row[j];
     358              if (jRow<=iRow+offset) {
     359                if (!used[jRow]) {
     360                  used[jRow]=1;
     361                  which[number++]=jRow;
     362                }
     363              }
     364            }
     365          }
     366        }
     367      } else {
     368        for (CoinBigIndex k=startRow;k<endRow;k++) {
     369          int iColumn=column[k];
     370          if (!whichDense_||!whichDense_[iColumn]) {
     371            CoinBigIndex start=columnStart[iColumn];
     372            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     373            for (CoinBigIndex j=start;j<end;j++) {
     374              int jRow=row[j];
     375              if (jRow>=iRow+offset) {
     376                if (!used[jRow]) {
     377                  used[jRow]=1;
     378                  which[number++]=jRow;
     379                }
     380              }
     381            }
     382          }
     383        }
     384      }
     385      sizeFactor_ += number;
     386      int j;
     387      for (j=0;j<number;j++)
     388        used[which[j]]=0;
     389    }
     390    delete [] which;
     391    // Now we have size - create arrays and fill in
     392    try {
     393      choleskyRow_ = new int [sizeFactor_];
     394    }
     395    catch (...) {
     396      // no memory
     397      delete [] choleskyStart_;
     398      choleskyStart_=NULL;
     399      return -1;
     400    }
     401    sizeFactor_=0;
     402    which = choleskyRow_;
     403    for (iRow=0;iRow<numberRows_;iRow++) {
     404      int number=0;
     405      // make sure diagonal exists if includeDiagonal
     406      if (!offset) {
     407        which[0]=iRow;
     408        used[iRow]=1;
     409        number=1;
     410      }
     411      choleskyStart_[iRow]=sizeFactor_;
     412      CoinBigIndex startRow=rowStart[iRow];
     413      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
     414      if (lowerTriangular) {
     415        for (CoinBigIndex k=startRow;k<endRow;k++) {
     416          int iColumn=column[k];
     417          if (!whichDense_||!whichDense_[iColumn]) {
     418            CoinBigIndex start=columnStart[iColumn];
     419            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     420            for (CoinBigIndex j=start;j<end;j++) {
     421              int jRow=row[j];
     422              if (jRow<=iRow+offset) {
     423                if (!used[jRow]) {
     424                  used[jRow]=1;
     425                  which[number++]=jRow;
     426                }
     427              }
     428            }
     429          }
     430        }
     431      } else {
     432        for (CoinBigIndex k=startRow;k<endRow;k++) {
     433          int iColumn=column[k];
     434          if (!whichDense_||!whichDense_[iColumn]) {
     435            CoinBigIndex start=columnStart[iColumn];
     436            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     437            for (CoinBigIndex j=start;j<end;j++) {
     438              int jRow=row[j];
     439              if (jRow>=iRow+offset) {
     440                if (!used[jRow]) {
     441                  used[jRow]=1;
     442                  which[number++]=jRow;
     443                }
     444              }
     445            }
     446          }
     447        }
     448      }
     449      sizeFactor_ += number;
     450      int j;
     451      for (j=0;j<number;j++)
     452        used[which[j]]=0;
     453      // Sort
     454      std::sort(which,which+number);
     455      // move which on
     456      which += number;
     457    }
     458    choleskyStart_[numberRows_]=sizeFactor_;
     459    delete [] used;
     460    return 0;
     461  } else {
     462    int numberRowsModel = model_->numberRows();
     463    int numberColumns = model_->numberColumns();
     464    int numberTotal = numberColumns + numberRowsModel;
     465    numberRows_ = 2*numberRowsModel+numberColumns;
     466    rowsDropped_ = new char [numberRows_];
     467    memset(rowsDropped_,0,numberRows_);
     468    numberRowsDropped_=0;
     469    CoinPackedMatrix * quadratic = NULL;
     470    ClpQuadraticObjective * quadraticObj =
     471      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     472    if (quadraticObj)
     473      quadratic = quadraticObj->quadraticObjective();
     474    int numberElements = model_->clpMatrix()->getNumElements();
     475    numberElements = numberElements+2*numberRowsModel+numberTotal;
     476    if (quadratic)
     477      numberElements += quadratic->getNumElements();
     478    // Space for starts
     479    choleskyStart_ = new CoinBigIndex[numberRows_+1];
     480    const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     481    const int * columnLength = model_->clpMatrix()->getVectorLengths();
     482    const int * row = model_->clpMatrix()->getIndices();
     483    //const double * element = model_->clpMatrix()->getElements();
     484    // Now we have size - create arrays and fill in
     485    try {
     486      choleskyRow_ = new int [numberElements];
     487    }
     488    catch (...) {
     489      // no memory
     490      delete [] choleskyStart_;
     491      choleskyStart_=NULL;
     492      return -1;
     493    }
     494    int iRow,iColumn;
     495 
     496    sizeFactor_=0;
     497    // matrix
     498    if (lowerTriangular) {
     499      if (!quadratic) {
     500        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     501          choleskyStart_[iColumn]=sizeFactor_;
     502          choleskyRow_[sizeFactor_++]=iColumn;
     503          CoinBigIndex start=columnStart[iColumn];
     504          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     505          if (!includeDiagonal)
     506            start++;
     507          for (CoinBigIndex j=start;j<end;j++) {
     508            choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
     509          }
     510        }
     511      } else {
     512        // Quadratic
     513        const int * columnQuadratic = quadratic->getIndices();
     514        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     515        const int * columnQuadraticLength = quadratic->getVectorLengths();
     516        //const double * quadraticElement = quadratic->getElements();
     517        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     518          choleskyStart_[iColumn]=sizeFactor_;
     519          if (includeDiagonal)
     520            choleskyRow_[sizeFactor_++]=iColumn;
     521          for (CoinBigIndex j=columnQuadraticStart[iColumn];
     522               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     523            int jColumn = columnQuadratic[j];
     524            if (jColumn>iColumn)
     525              choleskyRow_[sizeFactor_++]=jColumn;
     526          }
     527          CoinBigIndex start=columnStart[iColumn];
     528          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     529          for (CoinBigIndex j=start;j<end;j++) {
     530            choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
     531          }
     532        }
     533      }
     534      // slacks
     535      for (;iColumn<numberTotal;iColumn++) {
     536        choleskyStart_[iColumn]=sizeFactor_;
     537        if (includeDiagonal)
     538          choleskyRow_[sizeFactor_++]=iColumn;
     539        choleskyRow_[sizeFactor_++]=iColumn-numberColumns+numberTotal;
     540      }
     541      // Transpose - nonzero diagonal (may regularize)
     542      for (iRow=0;iRow<numberRowsModel;iRow++) {
     543        choleskyStart_[iRow+numberTotal]=sizeFactor_;
     544        // diagonal
     545        if (includeDiagonal)
     546          choleskyRow_[sizeFactor_++]=iRow+numberTotal;
     547      }
     548      choleskyStart_[numberRows_]=sizeFactor_;
     549    } else {
     550      // transpose
     551      ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
     552      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     553      const int * rowLength = rowCopy->getVectorLengths();
     554      const int * column = rowCopy->getIndices();
     555      if (!quadratic) {
     556        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     557          choleskyStart_[iColumn]=sizeFactor_;
     558          if (includeDiagonal)
     559            choleskyRow_[sizeFactor_++]=iColumn;
     560        }
     561      } else {
     562        // Quadratic
     563        // transpose
     564        CoinPackedMatrix quadraticT;
     565        quadraticT.reverseOrderedCopyOf(*quadratic);
     566        const int * columnQuadratic = quadraticT.getIndices();
     567        const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
     568        const int * columnQuadraticLength = quadraticT.getVectorLengths();
     569        //const double * quadraticElement = quadraticT.getElements();
     570        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     571          choleskyStart_[iColumn]=sizeFactor_;
     572          for (CoinBigIndex j=columnQuadraticStart[iColumn];
     573               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     574            int jColumn = columnQuadratic[j];
     575            if (jColumn<iColumn)
     576              choleskyRow_[sizeFactor_++]=jColumn;
     577          }
     578          if (includeDiagonal)
     579            choleskyRow_[sizeFactor_++]=iColumn;
     580        }
     581      }
     582      int iRow;
     583      // slacks
     584      for (iRow=0;iRow<numberRowsModel;iRow++) {
     585        choleskyStart_[iRow+numberColumns]=sizeFactor_;
     586        if (includeDiagonal)
     587          choleskyRow_[sizeFactor_++]=iRow+numberColumns;
     588      }
     589      for (iRow=0;iRow<numberRowsModel;iRow++) {
     590        choleskyStart_[iRow+numberTotal]=sizeFactor_;
     591        CoinBigIndex start=rowStart[iRow];
     592        CoinBigIndex end=rowStart[iRow]+rowLength[iRow];
     593        for (CoinBigIndex j=start;j<end;j++) {
     594          choleskyRow_[sizeFactor_++]=column[j];
     595        }
     596        // slack
     597        choleskyRow_[sizeFactor_++]=numberColumns+iRow;
     598        if (includeDiagonal)
     599          choleskyRow_[sizeFactor_++]=iRow+numberTotal;
     600      }
     601      choleskyStart_[numberRows_]=sizeFactor_;
     602    }
     603  }
     604  return 0;
     605}
     606/* Orders rows and saves pointer to matrix.and model */
     607int
     608ClpCholeskyBase::order(ClpInterior * model)
     609{
     610  model_=model;
     611  int numberRowsModel = model_->numberRows();
     612  int numberColumns = model_->numberColumns();
     613  int numberTotal = numberColumns + numberRowsModel;
     614  CoinPackedMatrix * quadratic = NULL;
     615  ClpQuadraticObjective * quadraticObj =
     616    (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     617  if (quadraticObj)
     618    quadratic = quadraticObj->quadraticObjective();
     619  if (!doKKT_) {
     620    numberRows_ = model->numberRows();
     621  } else {
     622    numberRows_ = 2*numberRowsModel+numberColumns;
     623  }
     624  rowsDropped_ = new char [numberRows_];
     625  // use for counts
     626  memset(rowsDropped_,0,numberRows_);
     627  numberRowsDropped_=0;
     628  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     629  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     630  const int * columnLength = model_->clpMatrix()->getVectorLengths();
     631  const int * row = model_->clpMatrix()->getIndices();
     632  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
     633  const int * rowLength = rowCopy_->getVectorLengths();
     634  const int * column = rowCopy_->getIndices();
     635  // We need two arrays for counts
     636  int * which = new int [numberRows_];
     637  int * used = new int[numberRows_+1];
     638  CoinZeroN(used,numberRows_);
     639  int iRow;
     640  sizeFactor_=0;
     641  permute_ = new int[numberRows_];
     642  for (iRow=0;iRow<numberRows_;iRow++)
     643    permute_[iRow]=iRow;
     644  if (!doKKT_) {
     645    int numberDense=0;
     646    if (denseThreshold_>0) {
     647      delete [] whichDense_;
     648      delete [] denseColumn_;
     649      delete dense_;
     650      whichDense_ = new char[numberColumns];
     651      int iColumn;
     652      used[numberRows_]=0;
     653      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     654        int length = columnLength[iColumn];
     655        used[length] += 1;
     656      }
     657      int nLong=0;
     658      int stop = CoinMax(denseThreshold_/2,100);
     659      for (iRow=numberRows_;iRow>=stop;iRow--) {
     660        if (used[iRow])
     661          printf("%d columns are of length %d\n",used[iRow],iRow);
     662        nLong += used[iRow];
     663        if (nLong>50||nLong>(numberColumns>>2))
     664          break;
     665      }
     666      CoinZeroN(used,numberRows_);
     667      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     668        if (columnLength[iColumn]<denseThreshold_) {
     669          whichDense_[iColumn]=0;
     670        } else {
     671          whichDense_[iColumn]=1;
     672          numberDense++;
     673        }
     674      }
     675      if (!numberDense||numberDense>100) {
     676        // free
     677        delete [] whichDense_;
     678        whichDense_=NULL;
     679        denseColumn_=NULL;
     680        dense_=NULL;
     681      } else {
     682        // space for dense columns
     683        denseColumn_ = new double [numberDense*numberRows_];
     684        // dense cholesky
     685        dense_ = new ClpCholeskyDense();
     686        dense_->reserveSpace(NULL,numberDense);
     687        printf("Taking %d columns as dense\n",numberDense);
     688      }
     689    }
     690    /*
     691       Get row counts and size
     692    */
     693    for (iRow=0;iRow<numberRows_;iRow++) {
     694      int number=1;
     695      // make sure diagonal exists
     696      which[0]=iRow;
     697      used[iRow]=1;
     698      CoinBigIndex startRow=rowStart[iRow];
     699      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
     700      for (CoinBigIndex k=startRow;k<endRow;k++) {
     701        int iColumn=column[k];
     702        if (!whichDense_||!whichDense_[iColumn]) {
     703          CoinBigIndex start=columnStart[iColumn];
     704          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     705          for (CoinBigIndex j=start;j<end;j++) {
     706            int jRow=row[j];
     707            if (jRow<iRow) {
     708              if (!used[jRow]) {
     709                used[jRow]=1;
     710                which[number++]=jRow;
     711                rowsDropped_[jRow]++;
     712              }
     713            }
     714          }
     715        }
     716    }
     717      sizeFactor_ += number;
     718      rowsDropped_[iRow]+=number;
     719      int j;
     720      for (j=0;j<number;j++)
     721        used[which[j]]=0;
     722    }
     723    CoinSort_2(rowsDropped_,rowsDropped_+numberRows_,permute_);
     724  } else {
     725    // KKT
     726    int numberElements = model_->clpMatrix()->getNumElements();
     727    numberElements = numberElements+2*numberRowsModel+numberTotal;
     728    if (quadratic)
     729      numberElements += quadratic->getNumElements();
     730    // off diagonal
     731    numberElements -= numberRows_;
     732    sizeFactor_=numberElements;
     733    // If we sort we need to redo symbolic etc
     734  }
     735  delete [] which;
     736  delete [] used;
     737  permuteInverse_ = new int [numberRows_];
     738  memset(rowsDropped_,0,numberRows_);
     739  for (iRow=0;iRow<numberRows_;iRow++) {
     740    //permute_[iRow]=iRow; // force no permute
     741    //permute_[iRow]=numberRows_-1-iRow; // force odd permute
     742    //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
     743    permuteInverse_[permute_[iRow]]=iRow;
     744  }
     745  return 0;
     746}
     747static int yyyyyy=0;
     748/* Does Symbolic factorization given permutation.
     749   This is called immediately after order.  If user provides this then
     750   user must provide factorize and solve.  Otherwise the default factorization is used
     751   returns non-zero if not enough memory */
     752int
     753ClpCholeskyBase::symbolic()
     754{
     755  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     756  const int * columnLength = model_->clpMatrix()->getVectorLengths();
     757  const int * row = model_->clpMatrix()->getIndices();
     758  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
     759  const int * rowLength = rowCopy_->getVectorLengths();
     760  const int * column = rowCopy_->getIndices();
     761  int numberRowsModel = model_->numberRows();
     762  int numberColumns = model_->numberColumns();
     763  int numberTotal = numberColumns + numberRowsModel;
     764  CoinPackedMatrix * quadratic = NULL;
     765  ClpQuadraticObjective * quadraticObj =
     766    (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     767  if (quadraticObj)
     768    quadratic = quadraticObj->quadraticObjective();
     769  // We need an array for counts
     770  int * used = new int[numberRows_+1];
     771  CoinZeroN(used,numberRows_);
     772  int iRow;
    104773  int iColumn;
    105   int numberColumns = model_->numberColumns();
    106   int numberTotal = numberRows_+numberColumns;
    107   double * region1Save = new double[numberTotal];
    108   for (iColumn=0;iColumn<numberTotal;iColumn++) {
    109     region1[iColumn] *= diagonal[iColumn];
    110     region1Save[iColumn]=region1[iColumn];
    111   }
    112   multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0);
    113   model_->clpMatrix()->times(1.0,region1,region2);
    114   double maximumRHS = maximumAbsElement(region2,numberRows_);
    115   double scale=1.0;
    116   double unscale=1.0;
    117   if (maximumRHS>1.0e-30) {
    118     if (maximumRHS<=0.5) {
    119       double factor=2.0;
    120       while (maximumRHS<=0.5) {
    121         maximumRHS*=factor;
    122         scale*=factor;
    123       } /* endwhile */
    124     } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
    125       double factor=0.5;
    126       while (maximumRHS>=2.0) {
    127         maximumRHS*=factor;
    128         scale*=factor;
    129       } /* endwhile */
    130     }
    131     unscale=diagonalScaleFactor/scale;
     774  bool noMemory=false;
     775  CoinBigIndex * Astart = new CoinBigIndex[numberRows_+1];
     776  int * Arow=NULL;
     777  try {
     778    Arow = new int [sizeFactor_];
     779  }
     780  catch (...) {
     781    // no memory
     782    delete [] Astart;
     783    return -1;
     784  }
     785  choleskyStart_ = new int[numberRows_+1];
     786  link_ = new int[numberRows_];
     787  workInteger_ = new CoinBigIndex[numberRows_];
     788  indexStart_ = new CoinBigIndex[numberRows_];
     789  clique_ = new int[numberRows_];
     790  // Redo so permuted upper triangular
     791  sizeFactor_=0;
     792  int * which = Arow;
     793  if (!doKKT_) {
     794    for (iRow=0;iRow<numberRows_;iRow++) {
     795      int number=0;
     796      int iOriginalRow = permute_[iRow];
     797      Astart[iRow]=sizeFactor_;
     798      CoinBigIndex startRow=rowStart[iOriginalRow];
     799      CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
     800      for (CoinBigIndex k=startRow;k<endRow;k++) {
     801        int iColumn=column[k];
     802        if (!whichDense_||!whichDense_[iColumn]) {
     803          CoinBigIndex start=columnStart[iColumn];
     804          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     805          for (CoinBigIndex j=start;j<end;j++) {
     806            int jRow=row[j];
     807            int jNewRow = permuteInverse_[jRow];
     808            if (jNewRow<iRow) {
     809              if (!used[jNewRow]) {
     810                used[jNewRow]=1;
     811                which[number++]=jNewRow;
     812              }
     813            }
     814          }
     815        }
     816      }
     817      sizeFactor_ += number;
     818      int j;
     819      for (j=0;j<number;j++)
     820        used[which[j]]=0;
     821      // Sort
     822      std::sort(which,which+number);
     823      // move which on
     824      which += number;
     825    }
    132826  } else {
    133     //effectively zero
    134     scale=0.0;
    135     unscale=0.0;
     827    // KKT
     828    // transpose
     829    ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
     830    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     831    const int * rowLength = rowCopy->getVectorLengths();
     832    const int * column = rowCopy->getIndices();
     833    // temp
     834    bool permuted=false;
     835    for (iRow=0;iRow<numberRows_;iRow++) {
     836      if (permute_[iRow]!=iRow) {
     837        permuted=true;
     838        break;
     839      }
     840    }
     841    // but fake it
     842    if (yyyyyy==1) {
     843      for (iRow=0;iRow<numberRows_;iRow++) {
     844        //permute_[iRow]=iRow; // force no permute
     845        permute_[iRow]=numberRows_-1-iRow; // force odd permute
     846        //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
     847        permuteInverse_[permute_[iRow]]=iRow;
     848      }
     849      permuted=true;
     850    } else if (yyyyyy==2) {
     851      for (iRow=0;iRow<numberRows_;iRow++) {
     852        permute_[iRow]=iRow; // force no permute
     853        //permute_[iRow]=numberRows_-1-iRow; // force odd permute
     854        //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
     855        permuteInverse_[permute_[iRow]]=iRow;
     856      }
     857    }
     858    if (permuted) {
     859      // Need to permute - ugly
     860      if (!quadratic) {
     861        for (iRow=0;iRow<numberRows_;iRow++) {
     862          Astart[iRow]=sizeFactor_;
     863          int iOriginalRow = permute_[iRow];
     864          if (iOriginalRow<numberColumns) {
     865            // A may be upper triangular by mistake
     866            iColumn=iOriginalRow;
     867            CoinBigIndex start=columnStart[iColumn];
     868            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     869            for (CoinBigIndex j=start;j<end;j++) {
     870              int kRow = row[j]+numberTotal;
     871              kRow=permuteInverse_[kRow];
     872              if (kRow<iRow)
     873                Arow[sizeFactor_++]=kRow;
     874            }
     875          } else if (iOriginalRow<numberTotal) {
     876            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
     877            if (kRow<iRow)
     878              Arow[sizeFactor_++]=kRow;
     879          } else {
     880            int kRow = iOriginalRow-numberTotal;
     881            CoinBigIndex start=rowStart[kRow];
     882            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
     883            for (CoinBigIndex j=start;j<end;j++) {
     884              int jRow = column[j];
     885              int jNewRow = permuteInverse_[jRow];
     886              if (jNewRow<iRow)
     887                Arow[sizeFactor_++]=jNewRow;
     888            }
     889            // slack - should it be permute
     890            kRow = permuteInverse_[kRow+numberColumns];
     891            if (kRow<iRow)
     892              Arow[sizeFactor_++]=kRow;
     893          }
     894          // Sort
     895          std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
     896        }
     897      } else {
     898        // quadratic
     899        // transpose
     900        CoinPackedMatrix quadraticT;
     901        quadraticT.reverseOrderedCopyOf(*quadratic);
     902        const int * columnQuadratic = quadraticT.getIndices();
     903        const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
     904        const int * columnQuadraticLength = quadraticT.getVectorLengths();
     905        for (iRow=0;iRow<numberRows_;iRow++) {
     906          Astart[iRow]=sizeFactor_;
     907          int iOriginalRow = permute_[iRow];
     908          if (iOriginalRow<numberColumns) {
     909            // Quadratic bit
     910            CoinBigIndex j;
     911            for ( j=columnQuadraticStart[iOriginalRow];
     912                  j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
     913              int jColumn = columnQuadratic[j];
     914              int jNewColumn = permuteInverse_[jColumn];
     915              if (jNewColumn<iRow)
     916                Arow[sizeFactor_++]=jNewColumn;
     917            }
     918            // A may be upper triangular by mistake
     919            iColumn=iOriginalRow;
     920            CoinBigIndex start=columnStart[iColumn];
     921            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     922            for (j=start;j<end;j++) {
     923              int kRow = row[j]+numberTotal;
     924              kRow=permuteInverse_[kRow];
     925              if (kRow<iRow)
     926                Arow[sizeFactor_++]=kRow;
     927            }
     928          } else if (iOriginalRow<numberTotal) {
     929            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
     930            if (kRow<iRow)
     931              Arow[sizeFactor_++]=kRow;
     932          } else {
     933            int kRow = iOriginalRow-numberTotal;
     934            CoinBigIndex start=rowStart[kRow];
     935            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
     936            for (CoinBigIndex j=start;j<end;j++) {
     937              int jRow = column[j];
     938              int jNewRow = permuteInverse_[jRow];
     939              if (jNewRow<iRow)
     940                Arow[sizeFactor_++]=jNewRow;
     941            }
     942            // slack - should it be permute
     943            kRow = permuteInverse_[kRow+numberColumns];
     944            if (kRow<iRow)
     945              Arow[sizeFactor_++]=kRow;
     946          }
     947          // Sort
     948          std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
     949        }
     950      }
     951    } else {
     952      if (!quadratic) {
     953        for (iRow=0;iRow<numberRows_;iRow++) {
     954          Astart[iRow]=sizeFactor_;
     955        }
     956      } else {
     957        // Quadratic
     958        // transpose
     959        CoinPackedMatrix quadraticT;
     960        quadraticT.reverseOrderedCopyOf(*quadratic);
     961        const int * columnQuadratic = quadraticT.getIndices();
     962        const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
     963        const int * columnQuadraticLength = quadraticT.getVectorLengths();
     964        //const double * quadraticElement = quadraticT.getElements();
     965        for (iRow=0;iRow<numberColumns;iRow++) {
     966          int iOriginalRow = permute_[iRow];
     967          Astart[iRow]=sizeFactor_;
     968          for (CoinBigIndex j=columnQuadraticStart[iOriginalRow];
     969               j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
     970            int jColumn = columnQuadratic[j];
     971            int jNewColumn = permuteInverse_[jColumn];
     972            if (jNewColumn<iRow)
     973              Arow[sizeFactor_++]=jNewColumn;
     974          }
     975        }
     976      }
     977      int iRow;
     978      // slacks
     979      for (iRow=0;iRow<numberRowsModel;iRow++) {
     980        Astart[iRow+numberColumns]=sizeFactor_;
     981      }
     982      for (iRow=0;iRow<numberRowsModel;iRow++) {
     983        Astart[iRow+numberTotal]=sizeFactor_;
     984        CoinBigIndex start=rowStart[iRow];
     985        CoinBigIndex end=rowStart[iRow]+rowLength[iRow];
     986        for (CoinBigIndex j=start;j<end;j++) {
     987          Arow[sizeFactor_++]=column[j];
     988        }
     989        // slack
     990        Arow[sizeFactor_++]=numberColumns+iRow;
     991      }
     992    }
     993  }
     994  Astart[numberRows_]=sizeFactor_;
     995  firstDense_=numberRows_;
     996  symbolic1(Astart,Arow);
     997  // Now fill in indices
     998  try {
     999    // too big
     1000    choleskyRow_ = new int[sizeFactor_];
     1001  }
     1002  catch (...) {
     1003    // no memory
     1004    noMemory=true;
    1361005  }
    137   multiplyAdd(NULL,numberRows_,0.0,region2,scale);
    138   solve(region2);
    139   multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
    140   multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns,0.0);
    141   CoinZeroN(region1,numberColumns);
    142   model_->clpMatrix()->transposeTimes(1.0,region2,region1);
    143   for (iColumn=0;iColumn<numberTotal;iColumn++)
    144     region1[iColumn] = region1[iColumn]*diagonal[iColumn]-region1Save[iColumn];
    145   delete [] region1Save;
     1006  if (!noMemory) {
     1007    // Do lower triangular
     1008    sizeFactor_=0;
     1009    int * which = Arow;
     1010    if (!doKKT_) {
     1011      for (iRow=0;iRow<numberRows_;iRow++) {
     1012        int number=0;
     1013        int iOriginalRow = permute_[iRow];
     1014        Astart[iRow]=sizeFactor_;
     1015        if (!rowsDropped_[iOriginalRow]) {
     1016          CoinBigIndex startRow=rowStart[iOriginalRow];
     1017          CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
     1018          for (CoinBigIndex k=startRow;k<endRow;k++) {
     1019            int iColumn=column[k];
     1020            if (!whichDense_||!whichDense_[iColumn]) {
     1021              CoinBigIndex start=columnStart[iColumn];
     1022              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1023              for (CoinBigIndex j=start;j<end;j++) {
     1024                int jRow=row[j];
     1025                int jNewRow = permuteInverse_[jRow];
     1026                if (jNewRow>iRow&&!rowsDropped_[jRow]) {
     1027                  if (!used[jNewRow]) {
     1028                    used[jNewRow]=1;
     1029                    which[number++]=jNewRow;
     1030                  }
     1031                }
     1032              }
     1033            }
     1034          }
     1035          sizeFactor_ += number;
     1036          int j;
     1037          for (j=0;j<number;j++)
     1038            used[which[j]]=0;
     1039          // Sort
     1040          std::sort(which,which+number);
     1041          // move which on
     1042          which += number;
     1043        }
     1044      }
     1045    } else {
     1046      // KKT
     1047      // temp
     1048      bool permuted=false;
     1049      for (iRow=0;iRow<numberRows_;iRow++) {
     1050        if (permute_[iRow]!=iRow) {
     1051          permuted=true;
     1052          break;
     1053        }
     1054      }
     1055      // but fake it
     1056      for (iRow=0;iRow<numberRows_;iRow++) {
     1057        //permute_[iRow]=iRow; // force no permute
     1058        //permuteInverse_[permute_[iRow]]=iRow;
     1059      }
     1060      if (permuted) {
     1061        // Need to permute - ugly
     1062        if (!quadratic) {
     1063          for (iRow=0;iRow<numberRows_;iRow++) {
     1064            Astart[iRow]=sizeFactor_;
     1065            int iOriginalRow = permute_[iRow];
     1066            if (iOriginalRow<numberColumns) {
     1067              iColumn=iOriginalRow;
     1068              CoinBigIndex start=columnStart[iColumn];
     1069              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1070              for (CoinBigIndex j=start;j<end;j++) {
     1071                int kRow = row[j]+numberTotal;
     1072                kRow=permuteInverse_[kRow];
     1073                if (kRow>iRow)
     1074                  Arow[sizeFactor_++]=kRow;
     1075              }
     1076            } else if (iOriginalRow<numberTotal) {
     1077              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
     1078              if (kRow>iRow)
     1079                Arow[sizeFactor_++]=kRow;
     1080            } else {
     1081              int kRow = iOriginalRow-numberTotal;
     1082              CoinBigIndex start=rowStart[kRow];
     1083              CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
     1084              for (CoinBigIndex j=start;j<end;j++) {
     1085                int jRow = column[j];
     1086                int jNewRow = permuteInverse_[jRow];
     1087                if (jNewRow>iRow)
     1088                  Arow[sizeFactor_++]=jNewRow;
     1089              }
     1090              // slack - should it be permute
     1091              kRow = permuteInverse_[kRow+numberColumns];
     1092              if (kRow>iRow)
     1093                Arow[sizeFactor_++]=kRow;
     1094            }
     1095            // Sort
     1096            std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
     1097          }
     1098        } else {
     1099          // quadratic
     1100          const int * columnQuadratic = quadratic->getIndices();
     1101          const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1102          const int * columnQuadraticLength = quadratic->getVectorLengths();
     1103          for (iRow=0;iRow<numberRows_;iRow++) {
     1104            Astart[iRow]=sizeFactor_;
     1105            int iOriginalRow = permute_[iRow];
     1106            if (iOriginalRow<numberColumns) {
     1107              // Quadratic bit
     1108              CoinBigIndex j;
     1109              for ( j=columnQuadraticStart[iOriginalRow];
     1110                    j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
     1111                int jColumn = columnQuadratic[j];
     1112                int jNewColumn = permuteInverse_[jColumn];
     1113                if (jNewColumn>iRow)
     1114                  Arow[sizeFactor_++]=jNewColumn;
     1115              }
     1116              iColumn=iOriginalRow;
     1117              CoinBigIndex start=columnStart[iColumn];
     1118              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1119              for (j=start;j<end;j++) {
     1120                int kRow = row[j]+numberTotal;
     1121                kRow=permuteInverse_[kRow];
     1122                if (kRow>iRow)
     1123                  Arow[sizeFactor_++]=kRow;
     1124              }
     1125            } else if (iOriginalRow<numberTotal) {
     1126              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
     1127              if (kRow>iRow)
     1128                Arow[sizeFactor_++]=kRow;
     1129            } else {
     1130              int kRow = iOriginalRow-numberTotal;
     1131              CoinBigIndex start=rowStart[kRow];
     1132              CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
     1133              for (CoinBigIndex j=start;j<end;j++) {
     1134                int jRow = column[j];
     1135                int jNewRow = permuteInverse_[jRow];
     1136                if (jNewRow>iRow)
     1137                  Arow[sizeFactor_++]=jNewRow;
     1138              }
     1139              // slack - should it be permute
     1140              kRow = permuteInverse_[kRow+numberColumns];
     1141              if (kRow>iRow)
     1142                Arow[sizeFactor_++]=kRow;
     1143            }
     1144            // Sort
     1145            std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
     1146          }
     1147        }
     1148      } else {
     1149        if (!quadratic) {
     1150          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1151            Astart[iColumn]=sizeFactor_;
     1152            CoinBigIndex start=columnStart[iColumn];
     1153            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1154            for (CoinBigIndex j=start;j<end;j++) {
     1155              Arow[sizeFactor_++]=row[j]+numberTotal;
     1156            }
     1157          }
     1158        } else {
     1159          // Quadratic
     1160          const int * columnQuadratic = quadratic->getIndices();
     1161          const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1162          const int * columnQuadraticLength = quadratic->getVectorLengths();
     1163          //const double * quadraticElement = quadratic->getElements();
     1164          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1165            Astart[iColumn]=sizeFactor_;
     1166            for (CoinBigIndex j=columnQuadraticStart[iColumn];
     1167                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1168              int jColumn = columnQuadratic[j];
     1169              if (jColumn>iColumn)
     1170                Arow[sizeFactor_++]=jColumn;
     1171            }
     1172            CoinBigIndex start=columnStart[iColumn];
     1173            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1174            for (CoinBigIndex j=start;j<end;j++) {
     1175              Arow[sizeFactor_++]=row[j]+numberTotal;
     1176            }
     1177          }
     1178        }
     1179        // slacks
     1180        for (iRow=0;iRow<numberRowsModel;iRow++) {
     1181          Astart[iRow+numberColumns]=sizeFactor_;
     1182          Arow[sizeFactor_++]=iRow+numberTotal;
     1183        }
     1184        // Transpose - nonzero diagonal (may regularize)
     1185        for (iRow=0;iRow<numberRowsModel;iRow++) {
     1186          Astart[iRow+numberTotal]=sizeFactor_;
     1187        }
     1188      }
     1189      Astart[numberRows_]=sizeFactor_;
     1190    }
     1191    symbolic2(Astart,Arow);
     1192    if (sizeIndex_<sizeFactor_) {
     1193      int * indices=NULL;
     1194      try {
     1195        indices = new int[sizeIndex_];
     1196      }
     1197      catch (...) {
     1198        // no memory
     1199        noMemory=true;
     1200      }
     1201      if (!noMemory)  {
     1202        memcpy(indices,choleskyRow_,sizeIndex_*sizeof(int));
     1203        delete [] choleskyRow_;
     1204        choleskyRow_=indices;
     1205      }
     1206    }
     1207  }
     1208  delete [] used;
     1209  // Use cholesky regions
     1210  delete [] Astart;
     1211  delete [] Arow;
     1212  double flops=0.0;
     1213  for (iRow=0;iRow<numberRows_;iRow++) {
     1214    int length = choleskyStart_[iRow+1]-choleskyStart_[iRow];
     1215    flops += ((double) length) * (length + 2.0);
     1216  }
     1217  if (model_->messageHandler()->logLevel()>0)
     1218    std::cout<<sizeFactor_<<" elements in sparse Cholesky, flop count "<<flops<<std::endl;
     1219  try {
     1220    sparseFactor_ = new double [sizeFactor_];
     1221    workDouble_ = new double[numberRows_];
     1222    diagonal_ = new double[numberRows_];
     1223  }
     1224  catch (...) {
     1225    // no memory
     1226    noMemory=true;
     1227  }
     1228  if (noMemory) {
     1229    delete [] choleskyRow_;
     1230    choleskyRow_=NULL;
     1231    delete [] choleskyStart_;
     1232    choleskyStart_=NULL;
     1233    delete [] permuteInverse_;
     1234    permuteInverse_ = NULL;
     1235    delete [] permute_;
     1236    permute_ = NULL;
     1237    delete [] choleskyStart_;
     1238    choleskyStart_ = NULL;
     1239    delete [] indexStart_;
     1240    indexStart_ = NULL;
     1241    delete [] link_;
     1242    link_ = NULL;
     1243    delete [] workInteger_;
     1244    workInteger_ = NULL;
     1245    delete [] sparseFactor_;
     1246    sparseFactor_ = NULL;
     1247    delete [] workDouble_;
     1248    workDouble_ = NULL;
     1249    delete [] diagonal_;
     1250    diagonal_ = NULL;
     1251    delete [] clique_;
     1252    clique_ = NULL;
     1253    return -1;
     1254  }
     1255  return 0;
    1461256}
     1257static int xxxxxx=2;
     1258int
     1259ClpCholeskyBase::symbolic1(const CoinBigIndex * Astart, const int * Arow)
     1260{
     1261  if (xxxxxx!=2) {
     1262    int iRow;
     1263    int put=0;
     1264    choleskyStart_[0]=put;
     1265    for (iRow=0;iRow<numberRows_;iRow++) {
     1266      put += numberRows_-iRow-1;
     1267      choleskyStart_[iRow+1]=put;
     1268      clique_[iRow]=-1;
     1269    }
     1270    sizeFactor_ = choleskyStart_[numberRows_];
     1271    return sizeFactor_;
     1272  }
     1273  int * marked = (int *) workInteger_;
     1274  int iRow;
     1275  // may not need to do this here but makes debugging easier
     1276  for (iRow=0;iRow<numberRows_;iRow++) {
     1277    marked[iRow]=-1;
     1278    link_[iRow]=-1;
     1279    choleskyStart_[iRow]=0; // counts
     1280  }
     1281  for (iRow=0;iRow<numberRows_;iRow++) {
     1282    marked[iRow]=iRow;
     1283    for (CoinBigIndex j=Astart[iRow];j<Astart[iRow+1];j++) {
     1284      int kRow=Arow[j];
     1285      while (marked[kRow] != iRow ) {
     1286        if (link_[kRow] <0 )
     1287          link_[kRow]=iRow;
     1288        choleskyStart_[kRow]++;
     1289        marked[kRow]=iRow;
     1290        kRow = link_[kRow];
     1291      }
     1292    }
     1293  }
     1294  sizeFactor_=0;
     1295  for (iRow=0;iRow<numberRows_;iRow++) {
     1296    int number = choleskyStart_[iRow];
     1297    choleskyStart_[iRow]=sizeFactor_;
     1298    sizeFactor_ += number;
     1299  }
     1300  choleskyStart_[numberRows_]=sizeFactor_;
     1301  return sizeFactor_;;
     1302}
     1303void
     1304ClpCholeskyBase::symbolic2(const CoinBigIndex * Astart, const int * Arow)
     1305{
     1306  if (xxxxxx!=2) {
     1307    int iRow;
     1308    int put=0;
     1309    choleskyStart_[0]=put;
     1310    for (iRow=0;iRow<numberRows_;iRow++) {
     1311      put += numberRows_-iRow-1;
     1312      choleskyStart_[iRow+1]=put;
     1313      clique_[iRow]=-1;
     1314    }
     1315    memcpy(indexStart_,choleskyStart_,numberRows_*sizeof(int));
     1316    sizeFactor_ = choleskyStart_[numberRows_];
     1317    sizeIndex_ = sizeFactor_;
     1318    put=0;
     1319   
     1320    xxxxxx=2;
     1321    for (iRow=0;iRow<numberRows_;iRow++) {
     1322      for (int j=iRow+1;j<numberRows_;j++) {
     1323        choleskyRow_[put++]=j;
     1324      }
     1325    }
     1326    return;
     1327  }
     1328  int iRow;
     1329#if 0
     1330  int * marked = (int *) workInteger_;
     1331  int * first = new int[numberRows_];
     1332  int * index = new int[numberRows_];
     1333  sizeFactor_=0;
     1334  sizeIndex_=0;
     1335  choleskyStart_[0]=0;
     1336  // may not need to do this here but makes debugging easier
     1337  for (iRow=0;iRow<numberRows_;iRow++) {
     1338    marked[iRow]=-1;
     1339    link_[iRow]=-1;
     1340    first[iRow]=Astart[iRow];
     1341  }
     1342  for (iRow=0;iRow<numberRows_;iRow++) {
     1343    int kRow=iRow;
     1344    int number=0;
     1345    // skip diagonal
     1346    int n=Astart[iRow+1]-Astart[iRow]-1;
     1347    choleskyStart_[iRow+1]=sizeFactor_+n;
     1348    memcpy(choleskyRow_+sizeFactor_,Arow+Astart[iRow]+1,n*sizeof(int));
     1349    first[iRow]=sizeFactor_;
     1350    while (kRow>=0) {
     1351      int nextRow=link_[kRow];
     1352      CoinBigIndex k=first[kRow];
     1353      CoinBigIndex end=choleskyStart_[kRow+1];
     1354      if (k+1<end) {
     1355        // update first
     1356        if (kRow!=iRow) {
     1357          assert (iRow==choleskyRow_[k]);
     1358          k++;
     1359          first[kRow]=k;
     1360        }
     1361        int jRow = choleskyRow_[k];
     1362        link_[kRow]=link_[jRow];
     1363        link_[jRow]=kRow;
     1364        for (;k<end;k++) {
     1365          int jRow = choleskyRow_[k];
     1366          if (marked[jRow]<0) {
     1367            // not marked yet
     1368            // may want to compare indices to last (later)
     1369            marked[jRow]=1;
     1370            index[number++]=jRow;
     1371          }
     1372        }
     1373      }
     1374      kRow=nextRow;
     1375    }
     1376    CoinBigIndex put=sizeFactor_;
     1377    // could look for clique
     1378    for (int i=0;i<number;i++) {
     1379      int jRow=index[i];
     1380      assert (jRow!=iRow);
     1381      choleskyRow_[put++]=jRow;
     1382      marked[jRow]=-1;
     1383    }
     1384    sizeFactor_ = put;
     1385    choleskyStart_[iRow+1]=sizeFactor_;
     1386    indexStart_[iRow]=sizeIndex_;
     1387    sizeIndex_ =put;
     1388  }
     1389  delete [] first;
     1390  delete [] index;
     1391  for (iRow=0;iRow<numberRows_;iRow++)
     1392    clique_[iRow]=-1;
     1393  return;
     1394#endif
     1395  int * mergeLink = clique_;
     1396  int * marker = (int *) workInteger_;
     1397  //int iRow;
     1398  for (iRow=0;iRow<numberRows_;iRow++) {
     1399    marker[iRow]=-1;
     1400    mergeLink[iRow]=-1;
     1401    link_[iRow]=-1; // not needed but makes debugging easier
     1402  }
     1403  int start=0;
     1404  int end=0;
     1405  choleskyStart_[0]=0;
     1406   
     1407  for (iRow=0;iRow<numberRows_;iRow++) {
     1408    int nz=0;
     1409    int merge = mergeLink[iRow];
     1410    bool marked=false;
     1411    if (merge<0)
     1412      marker[iRow]=iRow;
     1413    else
     1414      marker[iRow]=merge;
     1415    start = end;
     1416    int startSub=start;
     1417    link_[iRow]=numberRows_;
     1418    for (CoinBigIndex j=Astart[iRow];j<Astart[iRow+1];j++) {
     1419      int kRow=Arow[j];
     1420      int k=iRow;
     1421      int linked = link_[iRow];
     1422      while (linked<=kRow) {
     1423        k=linked;
     1424        linked = link_[k];
     1425      }
     1426      nz++;
     1427      link_[k]=kRow;
     1428      link_[kRow]=linked;
     1429      if (marker[kRow] != marker[iRow])
     1430        marked=true;
     1431    }
     1432    bool reuse=false;
     1433    // Check if we can re-use indices
     1434    if (!marked && merge>=0&&mergeLink[merge]<0&&0) {
     1435      // can re-use all
     1436      indexStart_[iRow]=indexStart_[merge]+1;
     1437      nz=choleskyStart_[merge+1]-(choleskyStart_[merge]+1);
     1438      reuse=true;
     1439    } else {
     1440      // See if we can re-use any
     1441      int k=mergeLink[iRow];
     1442      int maxLength=0;
     1443      while (k>=0) {
     1444        int length = choleskyStart_[k+1]-(choleskyStart_[k]+1);
     1445        int start = indexStart_[k]+1;
     1446        int stop = start+length;
     1447        if (length>maxLength) {
     1448          maxLength = length;
     1449          startSub = start;
     1450        }
     1451        int linked = iRow;
     1452        for (CoinBigIndex j=start;j<stop;j++) {
     1453          int kRow=choleskyRow_[j];
     1454          int kk=linked;
     1455          linked = link_[kk];
     1456          while (linked<kRow) {
     1457            kk=linked;
     1458            linked = link_[kk];
     1459          }
     1460          if (linked!=kRow) {
     1461            nz++;
     1462            link_[kk]=kRow;
     1463            link_[kRow]=linked;
     1464            linked=kRow;
     1465          }
     1466        }
     1467        k=mergeLink[k];
     1468      }
     1469      if (nz== maxLength)
     1470        reuse=true; // can re-use
     1471    }
     1472    reuse=false; //temp
     1473    if (!reuse) {
     1474      end += nz;
     1475      startSub=start;
     1476      int kRow = iRow;
     1477      for (int j=start;j<end;j++) {
     1478        kRow=link_[kRow];
     1479        choleskyRow_[j]=kRow;
     1480        assert (kRow<numberRows_);
     1481        marker[kRow]=iRow;
     1482      }
     1483      marker[iRow]=iRow;
     1484    }
     1485    indexStart_[iRow]=startSub;
     1486    choleskyStart_[iRow+1]=choleskyStart_[iRow]+nz;
     1487    if (nz>1) {
     1488      int kRow = choleskyRow_[startSub];
     1489      mergeLink[iRow]=mergeLink[kRow];
     1490      mergeLink[kRow]=iRow;
     1491    }
     1492    // should not be needed
     1493    std::sort(choleskyRow_+choleskyStart_[iRow]
     1494              ,choleskyRow_+choleskyStart_[iRow+1]);
     1495#ifdef CLP_DEBUG
     1496    int last=-1;
     1497    for (CoinBigIndex j=indexStart_[k];j<startSub;j++) {
     1498      int kRow=choleskyRow_[j];
     1499      assert (kRow>last);
     1500      last=kRow;
     1501    }
     1502#endif   
     1503  }
     1504  sizeFactor_ = choleskyStart_[numberRows_];
     1505  sizeIndex_ = start;
     1506  // find dense segment here
     1507  // Clean up clique info
     1508  for (iRow=0;iRow<numberRows_;iRow++)
     1509    clique_[iRow]=-1;
     1510  // should not be needed
     1511  for (iRow=0;iRow<numberRows_;iRow++) {
     1512    std::sort(choleskyRow_+choleskyStart_[iRow],
     1513              choleskyRow_+choleskyStart_[iRow+1]);
     1514  }
     1515}
     1516/* Factorize - filling in rowsDropped and returning number dropped */
     1517int
     1518ClpCholeskyBase::factorize(const double * diagonal, int * rowsDropped)
     1519{
     1520  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     1521  const int * columnLength = model_->clpMatrix()->getVectorLengths();
     1522  const int * row = model_->clpMatrix()->getIndices();
     1523  const double * element = model_->clpMatrix()->getElements();
     1524  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
     1525  const int * rowLength = rowCopy_->getVectorLengths();
     1526  const int * column = rowCopy_->getIndices();
     1527  const double * elementByRow = rowCopy_->getElements();
     1528  int numberColumns=model_->clpMatrix()->getNumCols();
     1529  //perturbation
     1530  double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     1531  perturbation=perturbation*perturbation;
     1532  if (perturbation>1.0) {
     1533    //if (model_->model()->logLevel()&4)
     1534      std::cout <<"large perturbation "<<perturbation<<std::endl;
     1535    perturbation=sqrt(perturbation);;
     1536    perturbation=1.0;
     1537  }
     1538  int iRow;
     1539  int iColumn;
     1540  double * work = workDouble_;
     1541  CoinZeroN(work,numberRows_);
     1542  int newDropped=0;
     1543  double largest=1.0;
     1544  double smallest=COIN_DBL_MAX;
     1545  int numberDense=0;
     1546  if (!doKKT_) {
     1547    const double * diagonalSlack = diagonal + numberColumns;
     1548    if (dense_)
     1549      numberDense=dense_->numberRows();
     1550    if (whichDense_) {
     1551      double * denseDiagonal = dense_->diagonal();
     1552      double * dense = denseColumn_;
     1553      int iDense=0;
     1554      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     1555        if (whichDense_[iColumn]) {
     1556          denseDiagonal[iDense++]=1.0/diagonal[iColumn];
     1557          CoinZeroN(dense,numberRows_);
     1558          CoinBigIndex start=columnStart[iColumn];
     1559          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1560          for (CoinBigIndex j=start;j<end;j++) {
     1561            int jRow=row[j];
     1562            dense[jRow] = element[j];
     1563          }
     1564          dense += numberRows_;
     1565        }
     1566      }
     1567    }
     1568    double delta2 = model_->delta(); // add delta*delta to diagonal
     1569    delta2 *= delta2;
     1570    // largest in initial matrix
     1571    double largest2=1.0e-20;
     1572    for (iRow=0;iRow<numberRows_;iRow++) {
     1573      double * put = sparseFactor_+choleskyStart_[iRow];
     1574      int * which = choleskyRow_+indexStart_[iRow];
     1575      int iOriginalRow = permute_[iRow];
     1576      int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
     1577      if (!rowLength[iOriginalRow])
     1578        rowsDropped_[iOriginalRow]=1;
     1579      if (!rowsDropped_[iOriginalRow]) {
     1580        CoinBigIndex startRow=rowStart[iOriginalRow];
     1581        CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
     1582        work[iRow] = diagonalSlack[iOriginalRow]+delta2;
     1583        for (CoinBigIndex k=startRow;k<endRow;k++) {
     1584          int iColumn=column[k];
     1585          if (!whichDense_||!whichDense_[iColumn]) {
     1586            CoinBigIndex start=columnStart[iColumn];
     1587            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1588            double multiplier = diagonal[iColumn]*elementByRow[k];
     1589            for (CoinBigIndex j=start;j<end;j++) {
     1590              int jRow=row[j];
     1591              int jNewRow = permuteInverse_[jRow];
     1592              if (jNewRow>=iRow&&!rowsDropped_[jRow]) {
     1593                double value=element[j]*multiplier;
     1594                work[jNewRow] += value;
     1595              }
     1596            }
     1597          }
     1598        }
     1599        diagonal_[iRow]=work[iRow];
     1600        largest2 = CoinMax(largest2,fabs(work[iRow]));
     1601        work[iRow]=0.0;
     1602        int j;
     1603        for (j=0;j<number;j++) {
     1604          int jRow = which[j];
     1605          put[j]=work[jRow];
     1606          largest2 = CoinMax(largest2,fabs(work[jRow]));
     1607          work[jRow]=0.0;
     1608        }
     1609      } else {
     1610        // dropped
     1611        diagonal_[iRow]=1.0;
     1612        int j;
     1613        for (j=1;j<number;j++) {
     1614          put[j]=0.0;
     1615        }
     1616      }
     1617    }
     1618    //check sizes
     1619    largest2*=1.0e-20;
     1620    largest = CoinMin(largest2,1.0e-11);
     1621    int numberDroppedBefore=0;
     1622    for (iRow=0;iRow<numberRows_;iRow++) {
     1623      int dropped=rowsDropped_[iRow];
     1624      // Move to int array
     1625      rowsDropped[iRow]=dropped;
     1626      if (!dropped) {
     1627        double diagonal = diagonal_[iRow];
     1628        if (diagonal>largest2) {
     1629          diagonal_[iRow]=diagonal+perturbation;
     1630        } else {
     1631          diagonal_[iRow]=diagonal+perturbation;
     1632          rowsDropped[iRow]=2;
     1633          numberDroppedBefore++;
     1634        }
     1635      }
     1636    }
     1637    doubleParameters_[10]=CoinMax(1.0e-20,largest);
     1638    integerParameters_[20]=0;
     1639    doubleParameters_[3]=0.0;
     1640    doubleParameters_[4]=COIN_DBL_MAX;
     1641    integerParameters_[34]=0; // say all must be positive
     1642    factorizePart2(rowsDropped);
     1643    newDropped=integerParameters_[20]+numberDroppedBefore;
     1644    largest=doubleParameters_[3];
     1645    smallest=doubleParameters_[4];
     1646    if (model_->messageHandler()->logLevel()>1)
     1647      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
     1648    choleskyCondition_=largest/smallest;
     1649    if (whichDense_) {
     1650      for (int i=0;i<numberRows_;i++) {
     1651        assert (diagonal_[i]>=0.0);
     1652        diagonal_[i]=sqrt(diagonal_[i]);
     1653      }
     1654      // Update dense columns (just L)
     1655      // Zero out dropped rows
     1656      int i;
     1657      for (i=0;i<numberDense;i++) {
     1658        double * a = denseColumn_+i*numberRows_;
     1659        for (int j=0;j<numberRows_;j++) {
     1660          if (rowsDropped[j])
     1661            a[j]=0.0;
     1662        }
     1663        solve(a,1);
     1664      }
     1665      dense_->resetRowsDropped();
     1666      longDouble * denseBlob = dense_->aMatrix();
     1667      double * denseDiagonal = dense_->diagonal();
     1668      // Update dense matrix
     1669      for (i=0;i<numberDense;i++) {
     1670        const double * a = denseColumn_+i*numberRows_;
     1671        // do diagonal
     1672        double value = denseDiagonal[i];
     1673        const double * b = denseColumn_+i*numberRows_;
     1674        for (int k=0;k<numberRows_;k++)
     1675          value += a[k]*b[k];
     1676        denseDiagonal[i]=value;
     1677        for (int j=i+1;j<numberDense;j++) {
     1678          double value = 0.0;
     1679          const double * b = denseColumn_+j*numberRows_;
     1680          for (int k=0;k<numberRows_;k++)
     1681            value += a[k]*b[k];
     1682          *denseBlob=value;
     1683          denseBlob++;
     1684        }
     1685      }
     1686      // dense cholesky (? long double)
     1687      int * dropped = new int [numberDense];
     1688      dense_->factorizePart2(dropped);
     1689      delete [] dropped;
     1690    }
     1691    bool cleanCholesky;
     1692    if (model_->numberIterations()<2000)
     1693      cleanCholesky=true;
     1694    else
     1695      cleanCholesky=false;
     1696    if (cleanCholesky) {
     1697      //drop fresh makes some formADAT easier
     1698      if (newDropped||numberRowsDropped_) {
     1699        newDropped=0;
     1700        for (int i=0;i<numberRows_;i++) {
     1701          char dropped = rowsDropped[i];
     1702          rowsDropped_[i]=dropped;
     1703          if (dropped==2) {
     1704            //dropped this time
     1705            rowsDropped[newDropped++]=i;
     1706            rowsDropped_[i]=0;
     1707          }
     1708        }
     1709        numberRowsDropped_=newDropped;
     1710        newDropped=-(2+newDropped);
     1711      }
     1712    } else {
     1713      if (newDropped) {
     1714        newDropped=0;
     1715        for (int i=0;i<numberRows_;i++) {
     1716          char dropped = rowsDropped[i];
     1717          rowsDropped_[i]=dropped;
     1718          if (dropped==2) {
     1719            //dropped this time
     1720            rowsDropped[newDropped++]=i;
     1721            rowsDropped_[i]=1;
     1722          }
     1723        }
     1724      }
     1725      numberRowsDropped_+=newDropped;
     1726      if (numberRowsDropped_&&0) {
     1727        std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
     1728          numberRowsDropped_<<" dropped)";
     1729        if (newDropped) {
     1730          std::cout<<" ( "<<newDropped<<" dropped this time)";
     1731        }
     1732        std::cout<<std::endl;
     1733      }
     1734    }
     1735  } else {
     1736    //KKT
     1737    CoinPackedMatrix * quadratic = NULL;
     1738    ClpQuadraticObjective * quadraticObj =
     1739      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     1740    if (quadraticObj)
     1741      quadratic = quadraticObj->quadraticObjective();
     1742    // matrix
     1743    int numberRowsModel = model_->numberRows();
     1744    int numberColumns = model_->numberColumns();
     1745    int numberTotal = numberColumns + numberRowsModel;
     1746    // temp
     1747    bool permuted=false;
     1748    for (iRow=0;iRow<numberRows_;iRow++) {
     1749      if (permute_[iRow]!=iRow) {
     1750        permuted=true;
     1751        break;
     1752      }
     1753    }
     1754    // but fake it
     1755    for (iRow=0;iRow<numberRows_;iRow++) {
     1756      //permute_[iRow]=iRow; // force no permute
     1757      //permuteInverse_[permute_[iRow]]=iRow;
     1758    }
     1759    if (permuted) {
     1760      double delta2 = model_->delta(); // add delta*delta to bottom
     1761      delta2 *= delta2;
     1762      // Need to permute - ugly
     1763      if (!quadratic) {
     1764        for (iRow=0;iRow<numberRows_;iRow++) {
     1765          double * put = sparseFactor_+choleskyStart_[iRow];
     1766          int * which = choleskyRow_+indexStart_[iRow];
     1767          int iOriginalRow = permute_[iRow];
     1768          if (iOriginalRow<numberColumns) {
     1769            iColumn=iOriginalRow;
     1770            double value = diagonal[iColumn];
     1771            if (fabs(value)>1.0e-100) {
     1772              value = 1.0/value;
     1773              largest = CoinMax(largest,fabs(value));
     1774              diagonal_[iRow] = -value;
     1775              CoinBigIndex start=columnStart[iColumn];
     1776              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1777              for (CoinBigIndex j=start;j<end;j++) {
     1778                int kRow = row[j]+numberTotal;
     1779                kRow=permuteInverse_[kRow];
     1780                if (kRow>iRow) {
     1781                  work[kRow]=element[j];
     1782                  largest = CoinMax(largest,fabs(element[j]));
     1783                }
     1784              }
     1785            } else {
     1786              diagonal_[iRow] = -value;
     1787            }
     1788          } else if (iOriginalRow<numberTotal) {
     1789            double value = diagonal[iOriginalRow];
     1790            if (fabs(value)>1.0e-100) {
     1791              value = 1.0/value;
     1792              largest = CoinMax(largest,fabs(value));
     1793            } else {
     1794              value = 1.0e100;
     1795            }
     1796            diagonal_[iRow] = -value;
     1797            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
     1798            if (kRow>iRow)
     1799              work[kRow]=-1.0;
     1800          } else {
     1801            diagonal_[iRow]=delta2;
     1802            int kRow = iOriginalRow-numberTotal;
     1803            CoinBigIndex start=rowStart[kRow];
     1804            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
     1805            for (CoinBigIndex j=start;j<end;j++) {
     1806              int jRow = column[j];
     1807              int jNewRow = permuteInverse_[jRow];
     1808              if (jNewRow>iRow) {
     1809                work[jNewRow]=elementByRow[j];
     1810                largest = CoinMax(largest,fabs(elementByRow[j]));
     1811              }
     1812            }
     1813            // slack - should it be permute
     1814            kRow = permuteInverse_[kRow+numberColumns];
     1815            if (kRow>iRow)
     1816              work[kRow]=-1.0;
     1817          }
     1818          CoinBigIndex j;
     1819          int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
     1820          for (j=0;j<number;j++) {
     1821            int jRow = which[j];
     1822            put[j]=work[jRow];
     1823            work[jRow]=0.0;
     1824          }
     1825        }
     1826      } else {
     1827        // quadratic
     1828        const int * columnQuadratic = quadratic->getIndices();
     1829        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1830        const int * columnQuadraticLength = quadratic->getVectorLengths();
     1831        const double * quadraticElement = quadratic->getElements();
     1832        for (iRow=0;iRow<numberRows_;iRow++) {
     1833          double * put = sparseFactor_+choleskyStart_[iRow];
     1834          int * which = choleskyRow_+indexStart_[iRow];
     1835          int iOriginalRow = permute_[iRow];
     1836          if (iOriginalRow<numberColumns) {
     1837            CoinBigIndex j;
     1838            iColumn=iOriginalRow;
     1839            double value = diagonal[iColumn];
     1840            if (fabs(value)>1.0e-100) {
     1841              value = 1.0/value;
     1842              for (j=columnQuadraticStart[iColumn];
     1843                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1844                int jColumn = columnQuadratic[j];
     1845                int jNewColumn = permuteInverse_[jColumn];
     1846                if (jNewColumn>iRow) {
     1847                  work[jNewColumn]=-quadraticElement[j];
     1848                } else if (iColumn==jColumn) {
     1849                  value += quadraticElement[j];
     1850                }
     1851              }
     1852              largest = CoinMax(largest,fabs(value));
     1853              diagonal_[iRow] = -value;
     1854              CoinBigIndex start=columnStart[iColumn];
     1855              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1856              for (j=start;j<end;j++) {
     1857                int kRow = row[j]+numberTotal;
     1858                kRow=permuteInverse_[kRow];
     1859                if (kRow>iRow) {
     1860                  work[kRow]=element[j];
     1861                  largest = CoinMax(largest,fabs(element[j]));
     1862                }
     1863              }
     1864            } else {
     1865              diagonal_[iRow] = -value;
     1866            }
     1867          } else if (iOriginalRow<numberTotal) {
     1868            double value = diagonal[iOriginalRow];
     1869            if (fabs(value)>1.0e-100) {
     1870              value = 1.0/value;
     1871              largest = CoinMax(largest,fabs(value));
     1872            } else {
     1873              value = 1.0e100;
     1874            }
     1875            diagonal_[iRow] = -value;
     1876            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
     1877            if (kRow>iRow)
     1878              work[kRow]=-1.0;
     1879          } else {
     1880            diagonal_[iRow]=delta2;
     1881            int kRow = iOriginalRow-numberTotal;
     1882            CoinBigIndex start=rowStart[kRow];
     1883            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
     1884            for (CoinBigIndex j=start;j<end;j++) {
     1885              int jRow = column[j];
     1886              int jNewRow = permuteInverse_[jRow];
     1887              if (jNewRow>iRow) {
     1888                work[jNewRow]=elementByRow[j];
     1889                largest = CoinMax(largest,fabs(elementByRow[j]));
     1890              }
     1891            }
     1892            // slack - should it be permute
     1893            kRow = permuteInverse_[kRow+numberColumns];
     1894            if (kRow>iRow)
     1895              work[kRow]=-1.0;
     1896          }
     1897          CoinBigIndex j;
     1898          int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
     1899          for (j=0;j<number;j++) {
     1900            int jRow = which[j];
     1901            put[j]=work[jRow];
     1902            work[jRow]=0.0;
     1903          }
     1904          for (j=0;j<numberRows_;j++)
     1905            assert (!work[j]);
     1906        }
     1907      }
     1908    } else {
     1909      if (!quadratic) {
     1910        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1911          double * put = sparseFactor_+choleskyStart_[iColumn];
     1912          int * which = choleskyRow_+indexStart_[iColumn];
     1913          double value = diagonal[iColumn];
     1914          if (fabs(value)>1.0e-100) {
     1915            value = 1.0/value;
     1916            largest = CoinMax(largest,fabs(value));
     1917            diagonal_[iColumn] = -value;
     1918            CoinBigIndex start=columnStart[iColumn];
     1919            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1920            for (CoinBigIndex j=start;j<end;j++) {
     1921              //choleskyRow_[numberElements]=row[j]+numberTotal;
     1922              //sparseFactor_[numberElements++]=element[j];
     1923              work[row[j]+numberTotal]=element[j];
     1924              largest = CoinMax(largest,fabs(element[j]));
     1925            }
     1926          } else {
     1927            diagonal_[iColumn] = -value;
     1928          }
     1929          CoinBigIndex j;
     1930          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
     1931          for (j=0;j<number;j++) {
     1932            int jRow = which[j];
     1933            put[j]=work[jRow];
     1934            work[jRow]=0.0;
     1935          }
     1936        }
     1937      } else {
     1938        // Quadratic
     1939        const int * columnQuadratic = quadratic->getIndices();
     1940        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1941        const int * columnQuadraticLength = quadratic->getVectorLengths();
     1942        const double * quadraticElement = quadratic->getElements();
     1943        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1944          double * put = sparseFactor_+choleskyStart_[iColumn];
     1945          int * which = choleskyRow_+indexStart_[iColumn];
     1946          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
     1947          double value = diagonal[iColumn];
     1948          CoinBigIndex j;
     1949          if (fabs(value)>1.0e-100) {
     1950            value = 1.0/value;
     1951            for (j=columnQuadraticStart[iColumn];
     1952                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1953              int jColumn = columnQuadratic[j];
     1954              if (jColumn>iColumn) {
     1955                work[jColumn]=-quadraticElement[j];
     1956              } else if (iColumn==jColumn) {
     1957                value += quadraticElement[j];
     1958              }
     1959            }
     1960            largest = CoinMax(largest,fabs(value));
     1961            diagonal_[iColumn] = -value;
     1962            CoinBigIndex start=columnStart[iColumn];
     1963            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     1964            for (j=start;j<end;j++) {
     1965              work[row[j]+numberTotal]=element[j];
     1966              largest = CoinMax(largest,fabs(element[j]));
     1967            }
     1968            for (j=0;j<number;j++) {
     1969              int jRow = which[j];
     1970              put[j]=work[jRow];
     1971              work[jRow]=0.0;
     1972            }
     1973          } else {
     1974            value = 1.0e100;
     1975            diagonal_[iColumn] = -value;
     1976            for (j=0;j<number;j++) {
     1977              int jRow = which[j];
     1978              put[j]=work[jRow];
     1979            }
     1980          }
     1981        }
     1982      }
     1983      // slacks
     1984      for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
     1985        double * put = sparseFactor_+choleskyStart_[iColumn];
     1986        int * which = choleskyRow_+indexStart_[iColumn];
     1987        double value = diagonal[iColumn];
     1988        if (fabs(value)>1.0e-100) {
     1989          value = 1.0/value;
     1990          largest = CoinMax(largest,fabs(value));
     1991        } else {
     1992          value = 1.0e100;
     1993        }
     1994        diagonal_[iColumn] = -value;
     1995        work[iColumn-numberColumns+numberTotal]=-1.0;
     1996        CoinBigIndex j;
     1997        int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
     1998        for (j=0;j<number;j++) {
     1999          int jRow = which[j];
     2000          put[j]=work[jRow];
     2001          work[jRow]=0.0;
     2002        }
     2003      }
     2004      // Finish diagonal
     2005      double delta2 = model_->delta(); // add delta*delta to bottom
     2006      delta2 *= delta2;
     2007      for (iRow=0;iRow<numberRowsModel;iRow++) {
     2008        double * put = sparseFactor_+choleskyStart_[iRow+numberTotal];
     2009        diagonal_[iRow+numberTotal]=delta2;
     2010        CoinBigIndex j;
     2011        int number = choleskyStart_[iRow+numberTotal+1]-choleskyStart_[iRow+numberTotal];
     2012        for (j=0;j<number;j++) {
     2013          put[j]=0.0;
     2014        }
     2015      }
     2016    }
     2017    //check sizes
     2018    largest*=1.0e-20;
     2019    largest = CoinMin(largest,1.0e-11);
     2020    doubleParameters_[10]=CoinMax(1.0e-20,largest);
     2021    integerParameters_[20]=0;
     2022    doubleParameters_[3]=0.0;
     2023    doubleParameters_[4]=COIN_DBL_MAX;
     2024    // Set up LDL cutoff
     2025    integerParameters_[34]=numberTotal;
     2026    // KKT
     2027    int * rowsDropped2 = new int[numberRows_];
     2028    CoinZeroN(rowsDropped2,numberRows_);
     2029    factorizePart2(rowsDropped2);
     2030    newDropped=integerParameters_[20];
     2031    largest=doubleParameters_[3];
     2032    smallest=doubleParameters_[4];
     2033    if (model_->messageHandler()->logLevel()>1)
     2034      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
     2035    choleskyCondition_=largest/smallest;
     2036    // Should save adjustments in ..R_
     2037    int n1=0,n2=0;
     2038    double * primalR = model_->primalR();
     2039    double * dualR = model_->dualR();
     2040    for (iRow=0;iRow<numberTotal;iRow++) {
     2041      if (rowsDropped2[iRow]) {
     2042        n1++;
     2043        //printf("row region1 %d dropped\n",iRow);
     2044        //rowsDropped_[iRow]=1;
     2045        rowsDropped_[iRow]=0;
     2046        primalR[iRow]=doubleParameters_[20];
     2047      } else {
     2048        rowsDropped_[iRow]=0;
     2049        primalR[iRow]=0.0;
     2050      }
     2051    }
     2052    for (;iRow<numberRows_;iRow++) {
     2053      if (rowsDropped2[iRow]) {
     2054        n2++;
     2055        //printf("row region2 %d dropped\n",iRow);
     2056        //rowsDropped_[iRow]=1;
     2057        rowsDropped_[iRow]=0;
     2058        dualR[iRow-numberTotal]=doubleParameters_[34];
     2059      } else {
     2060        rowsDropped_[iRow]=0;
     2061        dualR[iRow-numberTotal]=0.0;
     2062      }
     2063    }
     2064  }
     2065  status_=0;
     2066  return newDropped;
     2067}
     2068static /* Subroutine */ int ekkagmyldlt3(double *a, int *lda, int *n,
     2069                                  double *v, double *din)
     2070{
     2071    /* System generated locals */
     2072    int a_dim1, i__1, i__2;
     2073
     2074    /* Local variables */
     2075    int i, j, k;
     2076    double y, x0;
     2077    double t00;
     2078    double ai0, aj0, rd0;
     2079    double yj0;
     2080    double ddmin, ddmax;
     2081    double dsmall;
     2082
     2083/*      A contains the original matrix and the result is returned in lvals. */
     2084/*      A is garbage on output.  Both A and lvals are lower triangular normal.
     2085 */
     2086/*     do j = 0, N - 1 */
     2087/*       do i = 0, j - 1 */
     2088/*         x = v(i) * A(j,i) */
     2089/*         do k = j, N - 1 */
     2090/*           A(k,j) = A(k,j) - A(k,i) * x */
     2091/*         end do */
     2092/*       end do */
     2093/*       x = A(j,j) */
     2094/*       v(j) = x */
     2095/*       y = one/x */
     2096/*       A(j,j) = y */
     2097/*       do i = j + 1, N - 1 */
     2098/*         A(i,j) = A(i,j) * y */
     2099/*       end do */
     2100/*       print *,'Matrix after scalar iter = ',j */
     2101/*      call prtmat(A,lda,N) */
     2102/*     end do */
     2103/*     return */
     2104    /* Parameter adjustments */
     2105    a_dim1 = *lda;
     2106
     2107    /* Function Body */
     2108    dsmall = 1.0e-20;
     2109    ddmax = 0.0;
     2110    ddmin = 1.0e30;
     2111    i__1 = *n;
     2112    for (j = 0; j < i__1; j ++) {
     2113
     2114/* process column j */
     2115
     2116      //t00 = a[j + j * a_dim1];
     2117      t00 = din[j];
     2118      assert (t00==din[j]);
     2119      //assert (t00==v[j]);
     2120      for (k = 0; k <j; ++k) {
     2121        y = v[k];
     2122        aj0 = a[j + k * a_dim1];
     2123        yj0 = aj0 * y;
     2124        t00 -= aj0 * yj0;
     2125      }
     2126      x0 = fabs(t00);
     2127      if (x0 > ddmax) {
     2128        ddmax = t00;
     2129      }
     2130      if (x0 < ddmin) {
     2131        ddmin = t00;
     2132      }
     2133      if (x0 <= dsmall) {
     2134        abort();
     2135      }
     2136      rd0 = 1. / t00;
     2137      v[j] = t00;
     2138      //a[j + j * a_dim1] = rd0;
     2139      din[j]=rd0;
     2140
     2141/* process the remainder of column j */
     2142
     2143      i__2 = *n;
     2144      for (i = j+1; i < i__2; i ++) {
     2145        t00 = a[i + j * a_dim1];
     2146        for (k = 0; k <j; ++k) {
     2147          y = v[k];
     2148          aj0 = a[j + k * a_dim1] * y;
     2149          ai0 = a[i + k * a_dim1];
     2150          t00 -= ai0 * aj0;
     2151        }
     2152        t00 *= rd0;
     2153        a[i + j * a_dim1] = t00;
     2154      }
     2155    }
     2156    return 0;
     2157
     2158} /* ekkagmyldlt3 */
     2159/* Factorize - filling in rowsDropped and returning number dropped
     2160   in integerParam.
     2161*/
     2162void
     2163ClpCholeskyBase::factorizePart2(int * rowsDropped)
     2164{
     2165  double & largest=doubleParameters_[3];
     2166  double & smallest=doubleParameters_[4];
     2167  int numberDropped=integerParameters_[20];
     2168  // probably done before
     2169  largest=0.0;
     2170  smallest=COIN_DBL_MAX;
     2171  numberDropped=0;
     2172  double dropValue = doubleParameters_[10];
     2173  int firstPositive=integerParameters_[34];
     2174  double * d = ClpCopyOfArray(diagonal_,numberRows_);
     2175  int iRow;
     2176  if (0) {
     2177    double aa[40000];
     2178    //double v[200];
     2179    assert (numberRows_<200);
     2180    memset(d,0,numberRows_*sizeof(double));
     2181    memset(aa,0,numberRows_*numberRows_*sizeof(double));
     2182    for (iRow=0;iRow<numberRows_;iRow++) {
     2183      int base = numberRows_*iRow;
     2184      aa[base+iRow]=diagonal_[iRow];
     2185      for (int k=choleskyStart_[iRow];k<choleskyStart_[iRow+1];k++) {
     2186        int jRow = choleskyRow_[k];
     2187        aa[base+jRow]=sparseFactor_[k];
     2188      }
     2189    }
     2190    ekkagmyldlt3(aa, &numberRows_,&numberRows_,d,diagonal_);
     2191    //memcpy(diagonal_,v,numberRows_*sizeof(double));
     2192    for (iRow=0;iRow<numberRows_;iRow++) {
     2193      int base = numberRows_*iRow;
     2194      for (int k=choleskyStart_[iRow];k<choleskyStart_[iRow+1];k++) {
     2195        int jRow = choleskyRow_[k];
     2196        sparseFactor_[k]=aa[base+jRow];
     2197      }
     2198    }
     2199    return;
     2200  }
     2201  // minimum size before clique done
     2202#define MINCLIQUE INT_MAX
     2203  double * work = workDouble_;
     2204  CoinBigIndex * first = workInteger_;
     2205 
     2206  for (iRow=0;iRow<numberRows_;iRow++) {
     2207    link_[iRow]=-1;
     2208    work[iRow]=0.0;
     2209    first[iRow]=choleskyStart_[iRow];
     2210  }
     2211
     2212  int lastClique=-1;
     2213  bool inClique=false;
     2214  bool newClique=false;
     2215  bool endClique=false;
     2216 
     2217  for (iRow=0;iRow<firstDense_;iRow++) {
     2218    endClique=false;
     2219    if (clique_[iRow]>=0) {
     2220      // this is in a clique
     2221      inClique=true;
     2222      if (clique_[iRow]>lastClique) {
     2223        // new Clique
     2224        newClique=true;
     2225        // If we have clique going then signal to do old one
     2226        endClique=(lastClique>=0);
     2227      } else {
     2228        // Still in clique
     2229        newClique=false;
     2230      }
     2231    } else {
     2232      // not in clique
     2233      inClique=false;
     2234      newClique=false;
     2235      // If we have clique going then signal to do old one
     2236      endClique=(lastClique>=0);
     2237    }
     2238    lastClique=clique_[iRow];
     2239    if (endClique) {
     2240      // We have just finished updating a clique - do block pivot and clean up
     2241      abort(); //ekkchpp coding
     2242    }
     2243    if (newClique) {
     2244      // initialize new clique
     2245      abort();
     2246    }
     2247    // for each column L[*,kRow] that affects L[*,iRow]
     2248    double diagonalValue=diagonal_[iRow];
     2249    int nextRow = link_[iRow];
     2250    int kRow=0;
     2251    while (1) {
     2252      kRow=nextRow;
     2253      if (kRow<0)
     2254        break; // out of loop
     2255      nextRow=link_[kRow];
     2256      // Modify by outer product of L[*,irow] by L[*,krow] from first
     2257      CoinBigIndex k=first[kRow];
     2258      CoinBigIndex end=choleskyStart_[kRow+1];
     2259      assert(k<end);
     2260      double a_ik=sparseFactor_[k++];
     2261      double value1=d[kRow]*a_ik;
     2262      // update first
     2263      first[kRow]=k;
     2264      diagonalValue -= value1*a_ik;
     2265      if (k<end) {
     2266        int offset = indexStart_[kRow]-choleskyStart_[kRow];
     2267        int jRow = choleskyRow_[k+offset];
     2268        link_[kRow]=link_[jRow];
     2269        link_[jRow]=kRow;
     2270        if (clique_[kRow]<MINCLIQUE) {
     2271          for (;k<end;k++) {
     2272            int jRow = choleskyRow_[k+offset];
     2273            work[jRow] -= sparseFactor_[k]*value1;
     2274          }
     2275        } else {
     2276          // Clique
     2277          abort();
     2278        }
     2279      }
     2280    }
     2281    // Now apply
     2282    if (inClique) {
     2283      // in clique
     2284      abort();
     2285    } else {
     2286      // not in clique
     2287      int originalRow = permute_[iRow];
     2288      if (originalRow<firstPositive) {
     2289        // must be negative
     2290        if (diagonalValue<=-dropValue) {
     2291          smallest = CoinMin(smallest,-diagonalValue);
     2292          largest = CoinMax(largest,-diagonalValue);
     2293          d[iRow]=diagonalValue;
     2294          diagonalValue = 1.0/diagonalValue;
     2295        } else {
     2296          rowsDropped[originalRow]=2;
     2297          d[iRow]=-1.0e100;
     2298          diagonalValue=0.0;
     2299          integerParameters_[20]++;
     2300        }
     2301      } else {
     2302        // must be positive
     2303        if (diagonalValue>=dropValue) {
     2304          smallest = CoinMin(smallest,diagonalValue);
     2305          largest = CoinMax(largest,diagonalValue);
     2306          d[iRow]=diagonalValue;
     2307          diagonalValue = 1.0/diagonalValue;
     2308        } else {
     2309          rowsDropped[originalRow]=2;
     2310          d[iRow]=1.0e100;
     2311          diagonalValue=0.0;
     2312          integerParameters_[20]++;
     2313        }
     2314      }
     2315      diagonal_[iRow]=diagonalValue;
     2316      int offset = indexStart_[iRow]-choleskyStart_[iRow];
     2317      CoinBigIndex start = choleskyStart_[iRow];
     2318      CoinBigIndex end = choleskyStart_[iRow+1];
     2319      assert (first[iRow]==start);
     2320      if (start<end) {
     2321        int nextRow = choleskyRow_[start+offset];
     2322        link_[iRow]=link_[nextRow];
     2323        link_[nextRow]=iRow;
     2324        for (CoinBigIndex j=start;j<end;j++) {
     2325          int jRow = choleskyRow_[j+offset];
     2326          double value = sparseFactor_[j]+work[jRow];
     2327          work[jRow]=0.0;
     2328          sparseFactor_[j]= diagonalValue*value;
     2329        }
     2330      }
     2331    }
     2332  }
     2333 
     2334  delete [] d;
     2335  return;
     2336}
     2337/* Uses factorization to solve. */
     2338void
     2339ClpCholeskyBase::solve (double * region)
     2340{
     2341  if (!whichDense_) {
     2342    solve(region,3);
     2343  } else {
     2344    // dense columns
     2345    int i;
     2346    solve(region,1);
     2347    // do change;
     2348    int numberDense = dense_->numberRows();
     2349    double * change = new double[numberDense];
     2350    for (i=0;i<numberDense;i++) {
     2351      const double * a = denseColumn_+i*numberRows_;
     2352      double value =0.0;
     2353      for (int iRow=0;iRow<numberRows_;iRow++)
     2354        value += a[iRow]*region[iRow];
     2355      change[i]=value;
     2356    }
     2357    // solve
     2358    dense_->solve(change);
     2359    for (i=0;i<numberDense;i++) {
     2360      const double * a = denseColumn_+i*numberRows_;
     2361      double value = change[i];
     2362      for (int iRow=0;iRow<numberRows_;iRow++)
     2363        region[iRow] -= value*a[iRow];
     2364    }
     2365    delete [] change;
     2366    // and finish off
     2367    solve(region,2);
     2368  }
     2369}
     2370/* solve - 1 just first half, 2 just second half - 3 both.
     2371   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
     2372*/
     2373void
     2374ClpCholeskyBase::solve(double * region, int type)
     2375{
     2376  int i;
     2377  CoinBigIndex j;
     2378  for (i=0;i<numberRows_;i++) {
     2379    int iRow = permute_[i];
     2380    workDouble_[i] = region[iRow];
     2381  }
     2382  switch (type) {
     2383  case 1:
     2384    for (i=0;i<numberRows_;i++) {
     2385      double value=workDouble_[i];
     2386      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2387        int iRow = choleskyRow_[j];
     2388        workDouble_[iRow] -= sparseFactor_[j]*value;
     2389      }
     2390    }
     2391    for (i=0;i<numberRows_;i++) {
     2392      int iRow = permute_[i];
     2393      region[iRow]=workDouble_[i]*diagonal_[i];
     2394    }
     2395    break;
     2396  case 2:
     2397    for (i=numberRows_-1;i>=0;i--) {
     2398      double value=workDouble_[i]*diagonal_[i];
     2399      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2400        int iRow = choleskyRow_[j];
     2401        value -= sparseFactor_[j]*workDouble_[iRow];
     2402      }
     2403      workDouble_[i]=value;
     2404      int iRow = permute_[i];
     2405      region[iRow]=value;
     2406    }
     2407    break;
     2408  case 3:
     2409    for (i=0;i<numberRows_;i++) {
     2410      double value=workDouble_[i];
     2411      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2412        int iRow = choleskyRow_[j];
     2413        workDouble_[iRow] -= sparseFactor_[j]*value;
     2414      }
     2415    }
     2416    for (i=numberRows_-1;i>=0;i--) {
     2417      double value=workDouble_[i]*diagonal_[i];
     2418      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
     2419        int iRow = choleskyRow_[j];
     2420        value -= sparseFactor_[j]*workDouble_[iRow];
     2421      }
     2422      workDouble_[i]=value;
     2423      int iRow = permute_[i];
     2424      region[iRow]=value;
     2425    }
     2426    break;
     2427  }
     2428}
     2429
  • trunk/ClpCholeskyDense.cpp

    r320 r414  
    1010#include "ClpCholeskyDense.hpp"
    1111#include "ClpMessage.hpp"
     12#include "ClpQuadraticObjective.hpp"
    1213
    1314//#############################################################################
     
    2021ClpCholeskyDense::ClpCholeskyDense ()
    2122  : ClpCholeskyBase(),
    22     work_(NULL),
    23     rowCopy_(NULL)
     23    borrowSpace_(false)
    2424{
    2525  type_=11;;
     
    3030//-------------------------------------------------------------------
    3131ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs)
    32 : ClpCholeskyBase(rhs)
    33 {
    34   type_=rhs.type_;
    35   work_ = ClpCopyOfArray(rhs.work_,numberRows_*numberRows_);
    36   rowCopy_ = rhs.rowCopy_->clone();
     32  : ClpCholeskyBase(rhs),
     33    borrowSpace_(rhs.borrowSpace_)
     34{
     35  assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
    3736}
    3837
     
    4342ClpCholeskyDense::~ClpCholeskyDense ()
    4443{
    45   delete [] work_;
    46   delete rowCopy_;
     44  if (borrowSpace_) {
     45    // set NULL
     46    sparseFactor_=NULL;
     47    workDouble_=NULL;
     48    diagonal_=NULL;
     49  }
    4750}
    4851
     
    5457{
    5558  if (this != &rhs) {
     59    assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
    5660    ClpCholeskyBase::operator=(rhs);
    57     delete [] work_;
    58     work_ = ClpCopyOfArray(rhs.work_,numberRows_*numberRows_);
    59     delete rowCopy_;
    60     rowCopy_ = rhs.rowCopy_->clone();
     61    borrowSpace_=rhs.borrowSpace_;
    6162  }
    6263  return *this;
     
    6970  return new ClpCholeskyDense(*this);
    7071}
     72// If not power of 2 then need to redo a bit
     73#define BLOCK 16
     74#define BLOCKSHIFT 4
     75
     76#define BLOCKSQ BLOCK*BLOCK
     77#define BLOCKSQSHIFT BLOCKSHIFT+BLOCKSHIFT
     78#define number_blocks(x) (((x)+BLOCK-1)>>BLOCKSHIFT)
     79#define number_rows(x) ((x)<<BLOCKSHIFT)
     80#define number_entries(x) ((x)<<BLOCKSQSHIFT)
    7181/* Gets space */
    7282int
    73 ClpCholeskyDense::reserveSpace(int numberRows)
     83ClpCholeskyDense::reserveSpace(const ClpCholeskyBase * factor, int numberRows)
    7484{
    7585  numberRows_ = numberRows;
    76   work_ = new double [numberRows_*numberRows_];
    77   rowsDropped_ = new char [numberRows_];
    78   memset(rowsDropped_,0,numberRows_);
     86  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     87  // allow one stripe extra
     88  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
     89  sizeFactor_=numberBlocks*BLOCKSQ;
     90#define CHOL_COMPARE
     91#ifdef CHOL_COMPARE 
     92  sizeFactor_ += 95000;
     93#endif
     94  if (!factor) {
     95    sparseFactor_ = new double [sizeFactor_];
     96    rowsDropped_ = new char [numberRows_];
     97    memset(rowsDropped_,0,numberRows_);
     98    workDouble_ = new double[numberRows_];
     99    diagonal_ = new double[numberRows_];
     100  } else {
     101    borrowSpace_=true;
     102    int numberFull = factor->numberRows();
     103    sparseFactor_ = factor->sparseFactor()+(factor->size()-sizeFactor_);
     104    workDouble_ = factor->workDouble() + (numberFull-numberRows_);
     105    diagonal_ = factor->diagonal() + (numberFull-numberRows_);
     106  }
    79107  numberRowsDropped_=0;
    80108  return 0;
    81109}
     110/* Returns space needed */
     111CoinBigIndex
     112ClpCholeskyDense::space( int numberRows) const
     113{
     114 int numberBlocks = (numberRows+BLOCK-1)>>BLOCKSHIFT;
     115  // allow one stripe extra
     116  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
     117  CoinBigIndex sizeFactor=numberBlocks*BLOCKSQ;
     118#ifdef CHOL_COMPARE 
     119  sizeFactor += 95000;
     120#endif
     121  return sizeFactor;
     122 }
    82123/* Orders rows and saves pointer to matrix.and model */
    83124int
    84125ClpCholeskyDense::order(ClpInterior * model)
    85126{
    86   reserveSpace(model->numberRows());
    87127  model_=model;
     128  int numberRows;
     129  int numberRowsModel = model_->numberRows();
     130  int numberColumns = model_->numberColumns();
     131  if (!doKKT_) {
     132    numberRows = numberRowsModel;
     133  } else {
     134    numberRows = 2*numberRowsModel+numberColumns;
     135  }
     136  reserveSpace(NULL,numberRows);
    88137  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     138  return 0;
     139}
     140/* Does Symbolic factorization given permutation.
     141   This is called immediately after order.  If user provides this then
     142   user must provide factorize and solve.  Otherwise the default factorization is used
     143   returns non-zero if not enough memory */
     144int
     145ClpCholeskyDense::symbolic()
     146{
    89147  return 0;
    90148}
     
    93151ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped)
    94152{
     153  assert (!borrowSpace_);
    95154  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    96155  const int * columnLength = model_->clpMatrix()->getVectorLengths();
     
    102161  const double * elementByRow = rowCopy_->getElements();
    103162  int numberColumns=model_->clpMatrix()->getNumCols();
    104   CoinZeroN(work_,numberRows_*numberRows_);
     163  CoinZeroN(sparseFactor_,sizeFactor_);
     164  //perturbation
     165  double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
     166  perturbation=perturbation*perturbation;
     167  if (perturbation>1.0) {
     168    //if (model_->model()->logLevel()&4)
     169      std::cout <<"large perturbation "<<perturbation<<std::endl;
     170    perturbation=sqrt(perturbation);;
     171    perturbation=1.0;
     172  }
    105173  int iRow;
    106   double * work = work_;
    107   const double * diagonalSlack = diagonal + numberColumns;
    108   for (iRow=0;iRow<numberRows_;iRow++) {
    109     if (!rowsDropped_[iRow]) {
    110       CoinBigIndex startRow=rowStart[iRow];
    111       CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
    112       work[iRow] = diagonalSlack[iRow];
    113       for (CoinBigIndex k=startRow;k<endRow;k++) {
    114         int iColumn=column[k];
    115         CoinBigIndex start=columnStart[iColumn];
    116         CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    117         double multiplier = diagonal[iColumn]*elementByRow[k];
    118         for (CoinBigIndex j=start;j<end;j++) {
    119           int jRow=row[j];
    120           if (jRow>=iRow&&!rowsDropped_[jRow]) {
    121             double value=element[j]*multiplier;
    122             work[jRow] += value;
     174  int newDropped=0;
     175  double largest=1.0;
     176  double smallest=COIN_DBL_MAX;
     177  double delta2 = model_->delta(); // add delta*delta to diagonal
     178  delta2 *= delta2;
     179  if (!doKKT_) {
     180    double * work = sparseFactor_;
     181    work--; // skip diagonal
     182    int addOffset=numberRows_-1;
     183    const double * diagonalSlack = diagonal + numberColumns;
     184    // largest in initial matrix
     185    double largest2=1.0e-20;
     186    for (iRow=0;iRow<numberRows_;iRow++) {
     187      if (!rowsDropped_[iRow]) {
     188        CoinBigIndex startRow=rowStart[iRow];
     189        CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
     190        double diagonalValue = diagonalSlack[iRow]+delta2;
     191        for (CoinBigIndex k=startRow;k<endRow;k++) {
     192          int iColumn=column[k];
     193          CoinBigIndex start=columnStart[iColumn];
     194          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     195          double multiplier = diagonal[iColumn]*elementByRow[k];
     196          for (CoinBigIndex j=start;j<end;j++) {
     197            int jRow=row[j];
     198            if (!rowsDropped_[jRow]) {
     199              if (jRow>iRow) {
     200                double value=element[j]*multiplier;
     201                work[jRow] += value;
     202              } else if (jRow==iRow) {
     203                double value=element[j]*multiplier;
     204                diagonalValue += value;
     205              }
     206            }
    123207          }
     208        }
     209        for (int j=iRow+1;j<numberRows_;j++)
     210          largest2 = CoinMax(largest2,fabs(work[j]));
     211        diagonal_[iRow]=diagonalValue;
     212        largest2 = CoinMax(largest2,fabs(diagonalValue));
     213      } else {
     214        // dropped
     215        diagonal_[iRow]=1.0;
     216      }
     217      addOffset--;
     218      work += addOffset;
     219    }
     220    //check sizes
     221    largest2*=1.0e-20;
     222    largest = CoinMin(largest2,1.0e-11);
     223    int numberDroppedBefore=0;
     224    for (iRow=0;iRow<numberRows_;iRow++) {
     225      int dropped=rowsDropped_[iRow];
     226      // Move to int array
     227      rowsDropped[iRow]=dropped;
     228      if (!dropped) {
     229        double diagonal = diagonal_[iRow];
     230        if (diagonal>largest2) {
     231          diagonal_[iRow]=diagonal+perturbation;
     232        } else {
     233          diagonal_[iRow]=diagonal+perturbation;
     234          rowsDropped[iRow]=2;
     235          numberDroppedBefore++;
     236        }
     237      }
     238    }
     239    doubleParameters_[10]=CoinMax(1.0e-20,largest);
     240    integerParameters_[20]=0;
     241    doubleParameters_[3]=0.0;
     242    doubleParameters_[4]=COIN_DBL_MAX;
     243    integerParameters_[34]=0; // say all must be positive
     244#ifdef CHOL_COMPARE 
     245    if (numberRows_<200)
     246      factorizePart3(rowsDropped);
     247#endif
     248    factorizePart2(rowsDropped);
     249    newDropped=integerParameters_[20]+numberDroppedBefore;
     250    largest=doubleParameters_[3];
     251    smallest=doubleParameters_[4];
     252    if (model_->messageHandler()->logLevel()>1)
     253      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
     254    choleskyCondition_=largest/smallest;
     255    //drop fresh makes some formADAT easier
     256    if (newDropped||numberRowsDropped_) {
     257      newDropped=0;
     258      for (int i=0;i<numberRows_;i++) {
     259        char dropped = rowsDropped[i];
     260        rowsDropped_[i]=dropped;
     261        if (dropped==2) {
     262          //dropped this time
     263          rowsDropped[newDropped++]=i;
     264          rowsDropped_[i]=0;
     265        }
     266      }
     267      numberRowsDropped_=newDropped;
     268      newDropped=-(2+newDropped);
     269    }
     270  } else {
     271    // KKT
     272    CoinPackedMatrix * quadratic = NULL;
     273    ClpQuadraticObjective * quadraticObj =
     274      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
     275    if (quadraticObj)
     276      quadratic = quadraticObj->quadraticObjective();
     277    // matrix
     278    int numberRowsModel = model_->numberRows();
     279    int numberColumns = model_->numberColumns();
     280    int numberTotal = numberColumns + numberRowsModel;
     281    double * work = sparseFactor_;
     282    work--; // skip diagonal
     283    int addOffset=numberRows_-1;
     284    int iColumn;
     285    if (!quadratic) {
     286      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     287        double value = diagonal[iColumn];
     288        if (fabs(value)>1.0e-100) {
     289          value = 1.0/value;
     290          largest = CoinMax(largest,fabs(value));
     291          diagonal_[iColumn] = -value;
     292          CoinBigIndex start=columnStart[iColumn];
     293          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     294          for (CoinBigIndex j=start;j<end;j++) {
     295            //choleskyRow_[numberElements]=row[j]+numberTotal;
     296            //sparseFactor_[numberElements++]=element[j];
     297            work[row[j]+numberTotal]=element[j];
     298            largest = CoinMax(largest,fabs(element[j]));
     299          }
     300        } else {
     301          diagonal_[iColumn] = -value;
    124302        }
    125       }
    126     }
    127     work += numberRows_;
    128   }
    129   return factorizePart2(rowsDropped);
    130 }
    131 //#define CLP_DEBUG
     303        addOffset--;
     304        work += addOffset;
     305      }
     306    } else {
     307      // Quadratic
     308      const int * columnQuadratic = quadratic->getIndices();
     309      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     310      const int * columnQuadraticLength = quadratic->getVectorLengths();
     311      const double * quadraticElement = quadratic->getElements();
     312      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     313        double value = diagonal[iColumn];
     314        CoinBigIndex j;
     315        if (fabs(value)>1.0e-100) {
     316          value = 1.0/value;
     317          for (j=columnQuadraticStart[iColumn];
     318               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     319            int jColumn = columnQuadratic[j];
     320            if (jColumn>iColumn) {
     321              work[jColumn]=-quadraticElement[j];
     322            } else if (iColumn==jColumn) {
     323              value += quadraticElement[j];
     324            }
     325          }
     326          largest = CoinMax(largest,fabs(value));
     327          diagonal_[iColumn] = -value;
     328          CoinBigIndex start=columnStart[iColumn];
     329          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
     330          for (j=start;j<end;j++) {
     331            work[row[j]+numberTotal]=element[j];
     332            largest = CoinMax(largest,fabs(element[j]));
     333          }
     334        } else {
     335          value = 1.0e100;
     336          diagonal_[iColumn] = -value;
     337        }
     338        addOffset--;
     339        work += addOffset;
     340      }
     341    }
     342    // slacks
     343    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
     344      double value = diagonal[iColumn];
     345      if (fabs(value)>1.0e-100) {
     346        value = 1.0/value;
     347        largest = CoinMax(largest,fabs(value));
     348      } else {
     349        value = 1.0e100;
     350      }
     351      diagonal_[iColumn] = -value;
     352      work[iColumn-numberColumns+numberTotal]=-1.0;
     353      addOffset--;
     354      work += addOffset;
     355    }
     356    // Finish diagonal
     357    for (iRow=0;iRow<numberRowsModel;iRow++) {
     358      diagonal_[iRow+numberTotal]=delta2;
     359    }
     360    //check sizes
     361    largest*=1.0e-20;
     362    largest = CoinMin(largest,1.0e-11);
     363    doubleParameters_[10]=CoinMax(1.0e-20,largest);
     364    integerParameters_[20]=0;
     365    doubleParameters_[3]=0.0;
     366    doubleParameters_[4]=COIN_DBL_MAX;
     367    // Set up LDL cutoff
     368    integerParameters_[34]=numberTotal;
     369    // KKT
     370    int * rowsDropped2 = new int[numberRows_];
     371    CoinZeroN(rowsDropped2,numberRows_);
     372#ifdef CHOL_COMPARE 
     373    if (numberRows_<200)
     374      factorizePart3(rowsDropped2);
     375#endif
     376    factorizePart2(rowsDropped2);
     377    newDropped=integerParameters_[20];
     378    largest=doubleParameters_[3];
     379    smallest=doubleParameters_[4];
     380    if (model_->messageHandler()->logLevel()>1)
     381      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
     382    choleskyCondition_=largest/smallest;
     383    // Should save adjustments in ..R_
     384    int n1=0,n2=0;
     385    double * primalR = model_->primalR();
     386    double * dualR = model_->dualR();
     387    for (iRow=0;iRow<numberTotal;iRow++) {
     388      if (rowsDropped2[iRow]) {
     389        n1++;
     390        //printf("row region1 %d dropped\n",iRow);
     391        //rowsDropped_[iRow]=1;
     392        rowsDropped_[iRow]=0;
     393        primalR[iRow]=doubleParameters_[20];
     394      } else {
     395        rowsDropped_[iRow]=0;
     396        primalR[iRow]=0.0;
     397      }
     398    }
     399    for (;iRow<numberRows_;iRow++) {
     400      if (rowsDropped2[iRow]) {
     401        n2++;
     402        //printf("row region2 %d dropped\n",iRow);
     403        //rowsDropped_[iRow]=1;
     404        rowsDropped_[iRow]=0;
     405        dualR[iRow-numberTotal]=doubleParameters_[34];
     406      } else {
     407        rowsDropped_[iRow]=0;
     408        dualR[iRow-numberTotal]=0.0;
     409      }
     410    }
     411  }
     412  return 0;
     413}
    132414/* Factorize - filling in rowsDropped and returning number dropped */
    133 int
    134 ClpCholeskyDense::factorizePart2( int * rowsDropped)
     415void
     416ClpCholeskyDense::factorizePart3( int * rowsDropped)
    135417{
    136418  int iColumn;
    137   double * work = work_;
    138 #ifdef CLP_DEBUG
    139   double * save = NULL;
    140   if (numberRows_<20) {
    141     save = new double [ numberRows_*numberRows_];
    142     memcpy(save,work,numberRows_*numberRows_*sizeof(double));
    143   }
    144 #endif
     419  double * xx = sparseFactor_;
     420  double * yy = diagonal_;
     421  diagonal_ = sparseFactor_ + 40000;
     422  sparseFactor_=diagonal_ + numberRows_;
     423  //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
     424  memcpy(sparseFactor_,xx,40000*sizeof(double));
     425  memcpy(diagonal_,yy,numberRows_*sizeof(double));
    145426  int numberDropped=0;
     427  double largest=0.0;
     428  double smallest=COIN_DBL_MAX;
     429  double dropValue = doubleParameters_[10];
     430  int firstPositive=integerParameters_[34];
     431  double * work = sparseFactor_;
     432  // Allow for triangular
     433  int addOffset=numberRows_-1;
     434  work--;
    146435  for (iColumn=0;iColumn<numberRows_;iColumn++) {
    147436    int iRow;
    148     if (!rowsDropped_[iColumn]) {
    149       double diagonalValue = work[iColumn];
    150       if (diagonalValue>pivotTolerance_) {
    151         diagonalValue = sqrt(diagonalValue);
    152         work[iColumn]=diagonalValue;
    153         for (iRow=iColumn+1;iRow<numberRows_;iRow++)
    154           work[iRow] /= diagonalValue;
    155         double * work2 = work;
    156         for (int jColumn=iColumn+1;jColumn<numberRows_;jColumn++) {
    157           work2 += numberRows_;
    158           double value = work[jColumn];
    159           for (iRow=jColumn;iRow<numberRows_;iRow++)
    160           work2[iRow] -= value*work[iRow];
     437    int addOffsetNow = numberRows_-1;;
     438    double * workNow=sparseFactor_-1+iColumn;
     439    double diagonalValue = diagonal_[iColumn];
     440    for (iRow=0;iRow<iColumn;iRow++) {
     441      double aj = *workNow;
     442      addOffsetNow--;
     443      workNow += addOffsetNow;
     444      double multiplier = workDouble_[iRow];
     445      diagonalValue -= aj*aj*multiplier;
     446    }
     447    bool dropColumn=false;
     448    if (iColumn<firstPositive) {
     449      // must be negative
     450      if (diagonalValue<=-dropValue) {
     451        smallest = CoinMin(smallest,-diagonalValue);
     452        largest = CoinMax(largest,-diagonalValue);
     453        workDouble_[iColumn]=diagonalValue;
     454        diagonalValue = 1.0/diagonalValue;
     455      } else {
     456        dropColumn=true;
     457        workDouble_[iColumn]=-1.0e100;
     458        diagonalValue=0.0;
     459        integerParameters_[20]++;
     460      }
     461    } else {
     462      // must be positive
     463      if (diagonalValue>=dropValue) {
     464        smallest = CoinMin(smallest,diagonalValue);
     465        largest = CoinMax(largest,diagonalValue);
     466        workDouble_[iColumn]=diagonalValue;
     467        diagonalValue = 1.0/diagonalValue;
     468      } else {
     469        dropColumn=true;
     470        workDouble_[iColumn]=1.0e100;
     471        diagonalValue=0.0;
     472        integerParameters_[20]++;
     473      }
     474    }
     475    if (!dropColumn) {
     476      diagonal_[iColumn]=diagonalValue;
     477      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
     478        double value = work[iRow];
     479        workNow = sparseFactor_-1;
     480        int addOffsetNow = numberRows_-1;;
     481        for (int jColumn=0;jColumn<iColumn;jColumn++) {
     482          double aj = workNow[iColumn];
     483          double multiplier = workDouble_[jColumn];
     484          double ai = workNow[iRow];
     485          addOffsetNow--;
     486          workNow += addOffsetNow;
     487          value -= aj*ai*multiplier;
    161488        }
     489        work[iRow] = value*diagonalValue;
     490      }
     491    } else {
     492      // drop column
     493      rowsDropped[iColumn]=2;
     494      numberDropped++;
     495      diagonal_[iColumn]=0.0;
     496      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
     497        work[iRow] = 0.0;
     498      }
     499    }
     500    addOffset--;
     501    work += addOffset;
     502  }
     503  doubleParameters_[3]=largest;
     504  doubleParameters_[4]=smallest;
     505  integerParameters_[20]=numberDropped;
     506  sparseFactor_=xx;
     507  diagonal_=yy;
     508}
     509//#define POS_DEBUG
     510#ifdef POS_DEBUG
     511static int counter=0;
     512int ClpCholeskyDense::bNumber(const double * array, int &iRow, int &iCol)
     513{
     514  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     515  double * a = sparseFactor_+BLOCKSQ*numberBlocks;
     516  int k = array-a;
     517  assert ((k%BLOCKSQ)==0);
     518  iCol=0;
     519  int kk=k>>BLOCKSQSHIFT;
     520  //printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);
     521  k=kk;
     522  while (k>=numberBlocks) {
     523    iCol++;
     524    k -= numberBlocks;
     525    numberBlocks--;
     526  }
     527  iRow = iCol+k;
     528  counter++;
     529  if(counter>100000)
     530    exit(77);
     531  return kk;
     532 }
     533#endif
     534/* Factorize - filling in rowsDropped and returning number dropped */
     535void
     536ClpCholeskyDense::factorizePart2( int * rowsDropped)
     537{
     538  int iColumn;
     539  //double * xx = sparseFactor_;
     540  //double * yy = diagonal_;
     541  //diagonal_ = sparseFactor_ + 40000;
     542  //sparseFactor_=diagonal_ + numberRows_;
     543  //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
     544  //memcpy(sparseFactor_,xx,40000*sizeof(double));
     545  //memcpy(diagonal_,yy,numberRows_*sizeof(double));
     546  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     547  // later align on boundary
     548  double * a = sparseFactor_+BLOCKSQ*numberBlocks;
     549  int n=numberRows_;
     550  int nRound = numberRows_&(~(BLOCK-1));
     551  // adjust if exact
     552  if (nRound==n)
     553    nRound -= BLOCK;
     554  int sizeLastBlock=n-nRound;
     555  int get = n*(n-1)/2; // ? as no diagonal
     556  int block = numberBlocks*(numberBlocks+1)/2;
     557  int ifOdd;
     558  int rowLast;
     559  if (sizeLastBlock!=BLOCK) {
     560    double * aa = &a[(block-1)*BLOCKSQ];
     561    rowLast=nRound-1;
     562    ifOdd=1;
     563    int put = BLOCKSQ;
     564    // do last separately
     565    put -= (BLOCK-sizeLastBlock)*(BLOCK+1);
     566    for (iColumn=numberRows_-1;iColumn>=nRound;iColumn--) {
     567      int put2 = put;
     568      put -= BLOCK;
     569      for (int iRow=numberRows_-1;iRow>iColumn;iRow--) {
     570        aa[--put2] = sparseFactor_[--get];
     571        assert (aa+put2>=sparseFactor_+get);
     572      }
     573      // save diagonal as well
     574      aa[--put2]=diagonal_[iColumn];
     575    }
     576    n=nRound;
     577    block--;
     578  } else {
     579    // exact fit
     580    rowLast=numberRows_-1;
     581    ifOdd=0;
     582  }
     583  // Now main loop
     584  int nBlock=0;
     585  for (;n>0;n-=BLOCK) {
     586    double * aa = &a[(block-1)*BLOCKSQ];
     587    double * aaLast=NULL;
     588    int put = BLOCKSQ;
     589    int putLast=0;
     590    // see if we have small block
     591    if (ifOdd) {
     592      aaLast = &a[(block-1)*BLOCKSQ];
     593      aa = aaLast-BLOCKSQ;
     594      putLast = BLOCKSQ-BLOCK+sizeLastBlock;
     595    }
     596    for (iColumn=n-1;iColumn>=n-BLOCK;iColumn--) {
     597      if (aaLast) {
     598        // last bit
     599        for (int iRow=numberRows_-1;iRow>rowLast;iRow--) {
     600          aaLast[--putLast] = sparseFactor_[--get];
     601          assert (aaLast+putLast>=sparseFactor_+get);
     602        }
     603        putLast -= BLOCK-sizeLastBlock;
     604      }
     605      double * aPut = aa;
     606      int j=rowLast;
     607      for (int jBlock=0;jBlock<=nBlock;jBlock++) {
     608        int put2 = put;
     609        int last = CoinMax(j-BLOCK,iColumn);
     610        for (int iRow=j;iRow>last;iRow--) {
     611          aPut[--put2] = sparseFactor_[--get];
     612          assert (aPut+put2>=sparseFactor_+get);
     613        }
     614        if (j-BLOCK<iColumn) {
     615          // save diagonal as well
     616          aPut[--put2]=diagonal_[iColumn];
     617        }
     618        j -= BLOCK;
     619        aPut -=BLOCKSQ;
     620      }
     621      put -= BLOCK;
     622    }
     623    nBlock++;
     624    block -= nBlock+ifOdd;
     625  }
     626  factor(a,numberRows_,numberBlocks,
     627         diagonal_,workDouble_,rowsDropped);
     628  //sparseFactor_=xx;
     629  //diagonal_=yy;
     630}
     631// Non leaf recursive factor
     632void
     633ClpCholeskyDense::factor(double * a, int n, int numberBlocks,
     634            double * diagonal, double * work, int * rowsDropped)
     635{
     636  if (n<=BLOCK) {
     637    factorLeaf(a,n,diagonal,work,rowsDropped);
     638  } else {
     639    int nb=number_blocks((n+1)>>1);
     640    int nThis=number_rows(nb);
     641    double * aother;
     642    int nLeft=n-nThis;
     643    int nintri=(nb*(nb+1))>>1;
     644    int nbelow=(numberBlocks-nb)*nb;
     645    factor(a,nThis,numberBlocks,diagonal,work,rowsDropped);
     646    triRec(a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
     647    aother=a+number_entries(nintri+nbelow);
     648    recTri(a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
     649    factor(aother,nLeft,
     650        numberBlocks-nb,diagonal+nThis,work+nThis,rowsDropped);
     651  }
     652}
     653// Non leaf recursive triangle rectangle update
     654void
     655ClpCholeskyDense::triRec(double * aTri, int nThis, double * aUnder,
     656                         double * diagonal, double * work,
     657                         int nLeft, int iBlock, int jBlock,
     658                         int numberBlocks)
     659{
     660  if (nThis<=BLOCK&&nLeft<=BLOCK) {
     661    triRecLeaf(aTri,aUnder,diagonal,work,nLeft);
     662  } else if (nThis<nLeft) {
     663    int nb=number_blocks((nLeft+1)>>1);
     664    int nLeft2=number_rows(nb);
     665    triRec(aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
     666    triRec(aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
     667          iBlock+nb,jBlock,numberBlocks);
     668  } else {
     669    int nb=number_blocks((nThis+1)>>1);
     670    int nThis2=number_rows(nb);
     671    double * aother;
     672    int kBlock=jBlock+nb;
     673    int i;
     674    int nintri=(nb*(nb+1))>>1;
     675    int nbelow=(numberBlocks-nb)*nb;
     676    triRec(aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
     677    /* and rectangular update */
     678    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
     679      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
     680    aother=aUnder+number_entries(i);
     681    recRec(aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
     682          diagonal,work,iBlock,kBlock,jBlock,numberBlocks);
     683    triRec(aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
     684           work+nThis2,nLeft,
     685      iBlock-nb,kBlock-nb,numberBlocks-nb);
     686  }
     687}
     688// Non leaf recursive rectangle triangle update
     689void
     690ClpCholeskyDense::recTri(double * aUnder, int nTri, int nDo,
     691                         int iBlock, int jBlock,double * aTri,
     692                         double * diagonal, double * work,
     693                         int numberBlocks)
     694{
     695  if (nTri<=BLOCK&&nDo<=BLOCK) {
     696    recTriLeaf(aUnder,aTri,diagonal,work,nTri);
     697  } else if (nTri<nDo) {
     698    int nb=number_blocks((nDo+1)>>1);
     699    int nDo2=number_rows(nb);
     700    double * aother;
     701    int i;
     702    recTri(aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     703    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
     704      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
     705    aother=aUnder+number_entries(i);
     706    recTri(aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
     707           work+nDo2,numberBlocks-nb);
     708  } else {
     709    int nb=number_blocks((nTri+1)>>1);
     710    int nTri2=number_rows(nb);
     711    double * aother;
     712    int i;
     713    recTri(aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     714    /* and rectangular update */
     715    i=((numberBlocks-iBlock)*(numberBlocks-iBlock+1)-
     716      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb+1))>>1;
     717    aother=aTri+number_entries(nb);
     718    recRec(aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
     719           diagonal,work,iBlock+nb,iBlock,jBlock,numberBlocks);
     720    recTri(aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
     721         aTri+number_entries(i),diagonal,work,numberBlocks);
     722  }
     723}
     724/* Non leaf recursive rectangle rectangle update,
     725   nUnder is number of rows in iBlock,
     726   nUnderK is number of rows in kBlock
     727*/
     728void
     729ClpCholeskyDense::recRec(double * above, int nUnder, int nUnderK,
     730                         int nDo, double * aUnder, double *aOther, double * diagonal, double * work,
     731            int kBlock,int iBlock, int jBlock,
     732            int numberBlocks)
     733{
     734  if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) {
     735    recRecLeaf(above , aUnder ,  aOther, diagonal,work, nUnderK);
     736  } else if (nDo<=nUnderK&&nUnder<=nUnderK) {
     737    int nb=number_blocks((nUnderK+1)>>1);
     738    int nUnder2=number_rows(nb);
     739    recRec(above,nUnder,nUnder2,nDo,aUnder,aOther,diagonal,work,
     740               kBlock,iBlock,jBlock,numberBlocks);
     741    recRec(above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
     742        aOther+number_entries(nb),diagonal,work,kBlock+nb,iBlock,jBlock,numberBlocks);
     743  } else if (nUnderK<=nDo&&nUnder<=nDo) {
     744    int nb=number_blocks((nDo+1)>>1);
     745    int nDo2=number_rows(nb);
     746    int i;
     747    recRec(above,nUnder,nUnderK,nDo2,aUnder,aOther,diagonal,work,
     748         kBlock,iBlock,jBlock,numberBlocks);
     749    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
     750      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
     751    recRec(above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
     752           aUnder+number_entries(i),
     753           aOther,diagonal+nDo2,work+nDo2,
     754           kBlock-nb,iBlock-nb,jBlock,numberBlocks-nb);
     755  } else {
     756    int nb=number_blocks((nUnder+1)>>1);
     757    int nUnder2=number_rows(nb);
     758    int i;
     759    recRec(above,nUnder2,nUnderK,nDo,aUnder,aOther,diagonal,work,
     760               kBlock,iBlock,jBlock,numberBlocks);
     761    i=((numberBlocks-iBlock)*(numberBlocks-iBlock-1)-
     762      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb-1))>>1;
     763    recRec(above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
     764        aOther+number_entries(i),diagonal,work,kBlock,iBlock+nb,jBlock,numberBlocks);
     765  }
     766}
     767// Leaf recursive factor
     768void
     769ClpCholeskyDense::factorLeaf(double * a, int n,
     770                double * diagonal, double * work, int * rowsDropped)
     771{
     772  double largest=doubleParameters_[3];
     773  double smallest=doubleParameters_[4];
     774  double dropValue = doubleParameters_[10];
     775  int firstPositive=integerParameters_[34];
     776  int rowOffset=diagonal-diagonal_;
     777  int numberDropped=0;
     778  int i, j, k;
     779  double t00,temp1, * aa;
     780  aa = a-BLOCK;
     781  for (j = 0; j < n; j ++) {
     782    aa+=BLOCK;
     783    t00 = aa[j];
     784    for (k = 0; k < j; ++k) {
     785      double multiplier = work[k];
     786      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
     787    }
     788    bool dropColumn=false;
     789    double useT00=t00;
     790    if (j+rowOffset<firstPositive) {
     791      // must be negative
     792      if (t00<=-dropValue) {
     793        smallest = CoinMin(smallest,-t00);
     794        largest = CoinMax(largest,-t00);
     795        //aa[j]=t00;
     796        t00 = 1.0/t00;
    162797      } else {
    163         // drop column
    164         rowsDropped_[iColumn]=1;
    165         rowsDropped[numberDropped++]=iColumn;
    166         numberRowsDropped_++;
    167         // clean up as this is a debug version
    168         for (iRow=0;iRow<iColumn;iRow++)
    169           work_[iColumn+iRow*numberRows_]=0.0;
    170         for (iRow=iColumn;iRow<numberRows_;iRow++)
    171           work[iRow]=0.0;
     798        dropColumn=true;
     799        //aa[j]=-1.0e100;
     800        useT00=-1.0e-100;
     801        t00=0.0;
     802        integerParameters_[20]++;
    172803      }
    173804    } else {
    174       // clean up as this is a debug version
    175       for (iRow=0;iRow<iColumn;iRow++)
    176         work_[iColumn+iRow*numberRows_]=0.0;
    177       for (iRow=iColumn;iRow<numberRows_;iRow++)
    178         work[iRow]=0.0;
    179     }
    180     work += numberRows_; // move on pointer
    181   }
    182 #ifdef CLP_DEBUG
    183   if (save) {
    184     double * array = new double [numberRows_];
    185     int iColumn;
    186     for (iColumn=0;iColumn<numberRows_;iColumn++) {
    187       double * s = save;
    188       int iRow;
    189       for (iRow=0;iRow<iColumn;iRow++) {
    190         array[iRow]=s[iColumn];
    191         s += numberRows_;
    192       }
    193       for (iRow=iColumn;iRow<numberRows_;iRow++) {
    194         array[iRow]=s[iRow];
    195       }
    196       solve(array);
    197       for (iRow=0;iRow<numberRows_;iRow++) {
    198         if (iRow!=iColumn)
    199           assert (fabs(array[iRow])<1.0e-7);
    200         else
    201           assert (fabs(array[iRow]-1.0)<1.0e-7);
    202       }
    203     }
    204     delete [] array;
    205     delete [] save;
    206   }
    207 #endif
    208   return numberDropped;
     805      // must be positive
     806      if (t00>=dropValue) {
     807        smallest = CoinMin(smallest,t00);
     808        largest = CoinMax(largest,t00);
     809        //aa[j]=t00;
     810        t00 = 1.0/t00;
     811      } else {
     812        dropColumn=true;
     813        //aa[j]=1.0e100;
     814        useT00=1.0e-100;
     815        t00=0.0;
     816        integerParameters_[20]++;
     817      }
     818    }
     819    if (!dropColumn) {
     820      diagonal[j]=t00;
     821      work[j]=useT00;
     822      temp1 = t00;
     823      for (i=j+1;i<n;i++) {
     824        t00=aa[i];
     825        for (k = 0; k < j; ++k) {
     826          double multiplier = work[k];
     827          t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
     828        }
     829        aa[i]=t00*temp1;
     830      }
     831    } else {
     832      // drop column
     833      rowsDropped[j+rowOffset]=2;
     834      numberDropped++;
     835      diagonal[j]=0.0;
     836      //aa[j]=1.0e100;
     837      work[j]=1.0e100;
     838      for (i=j+1;i<n;i++) {
     839        aa[i]=0.0;
     840      }
     841    }
     842  }
     843  doubleParameters_[3]=largest;
     844  doubleParameters_[4]=smallest;
     845  integerParameters_[20]+=numberDropped;
     846}
     847// Leaf recursive triangle rectangle update
     848void
     849ClpCholeskyDense::triRecLeaf(double * aTri, double * aUnder, double * diagonal, double * work,
     850                int nUnder)
     851{
     852#ifdef POS_DEBUG
     853  int iru,icu;
     854  int iu=bNumber(aUnder,iru,icu);
     855  int irt,ict;
     856  int it=bNumber(aTri,irt,ict);
     857  //printf("%d %d\n",iu,it);
     858  printf("trirecleaf  under (%d,%d), tri (%d,%d)\n",
     859         iru,icu,irt,ict);
     860  assert (diagonal==diagonal_+ict*BLOCK);
     861#endif
     862  int j;
     863  double * aa;
     864#if BLOCK==16
     865  if (nUnder==BLOCK) {
     866    aa = aTri-2*BLOCK;
     867    for (j = 0; j < BLOCK; j +=2) {
     868      aa+=2*BLOCK;
     869      double temp0 = diagonal[j];
     870      double temp1 = diagonal[j+1];
     871      for (int i=0;i<BLOCK;i+=2) {
     872        double t00=aUnder[i+j*BLOCK];
     873        double t10=aUnder[i+ BLOCK + j*BLOCK];
     874        double t01=aUnder[i+1+j*BLOCK];
     875        double t11=aUnder[i+1+ BLOCK + j*BLOCK];
     876        for (int k = 0; k < j; ++k) {
     877          double multiplier=work[k];
     878          double au0 = aUnder[i + k * BLOCK]*multiplier;
     879          double au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
     880          double at0 = aTri[j + k * BLOCK];
     881          double at1 = aTri[j + 1 + k * BLOCK];
     882          t00 -= au0 * at0;
     883          t10 -= au0 * at1;
     884          t01 -= au1 * at0;
     885          t11 -= au1 * at1;
     886        }
     887        t00 *= temp0;
     888        double at1 = aTri[j + 1 + j*BLOCK]*work[j];
     889        t10 -= t00 * at1;
     890        t01 *= temp0;
     891        t11 -= t01 * at1;
     892        aUnder[i+j*BLOCK]=t00;
     893        aUnder[i+1+j*BLOCK]=t01;
     894        aUnder[i+ BLOCK + j*BLOCK]=t10*temp1;
     895        aUnder[i+1+ BLOCK + j*BLOCK]=t11*temp1;
     896      }
     897    }
     898  } else {
     899#endif
     900    aa = aTri-BLOCK;
     901    for (j = 0; j < BLOCK; j ++) {
     902      aa+=BLOCK;
     903      double temp1 = diagonal[j];
     904      for (int i=0;i<nUnder;i++) {
     905        double t00=aUnder[i+j*BLOCK];
     906        for (int k = 0; k < j; ++k) {
     907          double multiplier=work[k];
     908          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
     909        }
     910        aUnder[i+j*BLOCK]=t00*temp1;
     911      }
     912    }
     913#if BLOCK==16
     914  }
     915#endif
     916}
     917// Leaf recursive rectangle triangle update
     918void ClpCholeskyDense::recTriLeaf(double * aUnder, double * aTri,
     919                                  double * diagonal, double * work, int nUnder)
     920{
     921#ifdef POS_DEBUG
     922  int iru,icu;
     923  int iu=bNumber(aUnder,iru,icu);
     924  int irt,ict;
     925  int it=bNumber(aTri,irt,ict);
     926  //printf("%d %d\n",iu,it);
     927  printf("rectrileaf  under (%d,%d), tri (%d,%d)\n",
     928         iru,icu,irt,ict);
     929  assert (diagonal==diagonal_+icu*BLOCK);
     930#endif
     931  int i, j, k;
     932  double t00, * aa;
     933#if BLOCK==16
     934  if (nUnder==BLOCK) {
     935    double * aUnder2 = aUnder-2;
     936    aa = aTri-2*BLOCK;
     937    for (j = 0; j < BLOCK; j +=2) {
     938      double t00,t01,t10,t11;
     939      aa+=2*BLOCK;
     940      aUnder2+=2;
     941      t00=aa[j];
     942      t01=aa[j+1];
     943      t10=aa[j+1+BLOCK];
     944      for (k = 0; k < BLOCK; ++k) {
     945        double multiplier = work[k];
     946        double a0 = aUnder2[k * BLOCK];
     947        double a1 = aUnder2[1 + k * BLOCK];
     948        double x0 = a0 * multiplier;
     949        double x1 = a1 * multiplier;
     950        t00 -= a0 * x0;
     951        t01 -= a1 * x0;
     952        t10 -= a1 * x1;
     953      }
     954      aa[j]=t00;
     955      aa[j+1]=t01;
     956      aa[j+1+BLOCK]=t10;
     957      for (i=j+2;i< BLOCK;i+=2) {
     958        t00=aa[i];
     959        t01=aa[i+BLOCK];
     960        t10=aa[i+1];
     961        t11=aa[i+1+BLOCK];
     962        for (k = 0; k < BLOCK; ++k) {
     963          double multiplier = work[k];
     964          double a0 = aUnder2[k * BLOCK]*multiplier;
     965          double a1 = aUnder2[1 + k * BLOCK]*multiplier;
     966          t00 -= aUnder[i + k * BLOCK] * a0;
     967          t01 -= aUnder[i + k * BLOCK] * a1;
     968          t10 -= aUnder[i + 1 + k * BLOCK] * a0;
     969          t11 -= aUnder[i + 1 + k * BLOCK] * a1;
     970        }
     971        aa[i]=t00;
     972        aa[i+BLOCK]=t01;
     973        aa[i+1]=t10;
     974        aa[i+1+BLOCK]=t11;
     975      }
     976    }
     977  } else {
     978#endif
     979    aa = aTri-BLOCK;
     980    for (j = 0; j < nUnder; j ++) {
     981      aa+=BLOCK;
     982      for (i=j;i< nUnder;i++) {
     983        t00=aa[i];
     984        for (k = 0; k < BLOCK; ++k) {
     985          double multiplier = work[k];
     986          t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
     987        }
     988        aa[i]=t00;
     989      }
     990    }
     991#if BLOCK==16
     992  }
     993#endif
     994}
     995/* Leaf recursive rectangle rectangle update,
     996   nUnder is number of rows in iBlock,
     997   nUnderK is number of rows in kBlock
     998*/
     999void
     1000ClpCholeskyDense::recRecLeaf(double * above,
     1001                             double * aUnder, double *aOther, double * diagonal, double * work,
     1002                             int nUnder)
     1003{
     1004#ifdef POS_DEBUG
     1005  int ira,ica;
     1006  int ia=bNumber(above,ira,ica);
     1007  int iru,icu;
     1008  int iu=bNumber(aUnder,iru,icu);
     1009  int iro,ico;
     1010  int io=bNumber(aOther,iro,ico);
     1011  //printf("%d %d %d\n",ia,iu,io);
     1012  printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n",
     1013         ira,ica,iru,icu,iro,ico);
     1014  assert (diagonal==diagonal_+ica*BLOCK);
     1015#endif
     1016  int i, j, k;
     1017  double * aa;
     1018#if BLOCK==16
     1019  aa = aOther-4*BLOCK;
     1020  if (nUnder==BLOCK) {
     1021    for (j = 0; j < BLOCK; j +=4) {
     1022      aa+=4*BLOCK;
     1023      for (i=0;i< BLOCK;i+=4) {
     1024        double t00=aa[i+0+0*BLOCK];
     1025        double t10=aa[i+0+1*BLOCK];
     1026        double t20=aa[i+0+2*BLOCK];
     1027        double t30=aa[i+0+3*BLOCK];
     1028        double t01=aa[i+1+0*BLOCK];
     1029        double t11=aa[i+1+1*BLOCK];
     1030        double t21=aa[i+1+2*BLOCK];
     1031        double t31=aa[i+1+3*BLOCK];
     1032        double t02=aa[i+2+0*BLOCK];
     1033        double t12=aa[i+2+1*BLOCK];
     1034        double t22=aa[i+2+2*BLOCK];
     1035        double t32=aa[i+2+3*BLOCK];
     1036        double t03=aa[i+3+0*BLOCK];
     1037        double t13=aa[i+3+1*BLOCK];
     1038        double t23=aa[i+3+2*BLOCK];
     1039        double t33=aa[i+3+3*BLOCK];
     1040        for (k=0;k<BLOCK;k++) {
     1041          double multiplier = work[k];
     1042          double a00=aUnder[i+0+k*BLOCK]*multiplier;
     1043          double a01=aUnder[i+1+k*BLOCK]*multiplier;
     1044          double a02=aUnder[i+2+k*BLOCK]*multiplier;
     1045          double a03=aUnder[i+3+k*BLOCK]*multiplier;
     1046          t00 -= a00 * above[j + 0 + k * BLOCK];
     1047          t10 -= a00 * above[j + 1 + k * BLOCK];
     1048          t20 -= a00 * above[j + 2 + k * BLOCK];
     1049          t30 -= a00 * above[j + 3 + k * BLOCK];
     1050          t01 -= a01 * above[j + 0 + k * BLOCK];
     1051          t11 -= a01 * above[j + 1 + k * BLOCK];
     1052          t21 -= a01 * above[j + 2 + k * BLOCK];
     1053          t31 -= a01 * above[j + 3 + k * BLOCK];
     1054          t02 -= a02 * above[j + 0 + k * BLOCK];
     1055          t12 -= a02 * above[j + 1 + k * BLOCK];
     1056          t22 -= a02 * above[j + 2 + k * BLOCK];
     1057          t32 -= a02 * above[j + 3 + k * BLOCK];
     1058          t03 -= a03 * above[j + 0 + k * BLOCK];
     1059          t13 -= a03 * above[j + 1 + k * BLOCK];
     1060          t23 -= a03 * above[j + 2 + k * BLOCK];
     1061          t33 -= a03 * above[j + 3 + k * BLOCK];
     1062        }
     1063        aa[i+0+0*BLOCK]=t00;
     1064        aa[i+0+1*BLOCK]=t10;
     1065        aa[i+0+2*BLOCK]=t20;
     1066        aa[i+0+3*BLOCK]=t30;
     1067        aa[i+1+0*BLOCK]=t01;
     1068        aa[i+1+1*BLOCK]=t11;
     1069        aa[i+1+2*BLOCK]=t21;
     1070        aa[i+1+3*BLOCK]=t31;
     1071        aa[i+2+0*BLOCK]=t02;
     1072        aa[i+2+1*BLOCK]=t12;
     1073        aa[i+2+2*BLOCK]=t22;
     1074        aa[i+2+3*BLOCK]=t32;
     1075        aa[i+3+0*BLOCK]=t03;
     1076        aa[i+3+1*BLOCK]=t13;
     1077        aa[i+3+2*BLOCK]=t23;
     1078        aa[i+3+3*BLOCK]=t33;
     1079      }
     1080    }
     1081  } else {
     1082    int odd=nUnder&1;
     1083    int n=nUnder-odd;
     1084    for (j = 0; j < BLOCK; j +=4) {
     1085      aa+=4*BLOCK;
     1086      for (i=0;i< n;i+=2) {
     1087        double t00=aa[i+0*BLOCK];
     1088        double t10=aa[i+1*BLOCK];
     1089        double t20=aa[i+2*BLOCK];
     1090        double t30=aa[i+3*BLOCK];
     1091        double t01=aa[i+1+0*BLOCK];
     1092        double t11=aa[i+1+1*BLOCK];
     1093        double t21=aa[i+1+2*BLOCK];
     1094        double t31=aa[i+1+3*BLOCK];
     1095        for (k=0;k<BLOCK;k++) {
     1096          double multiplier = work[k];
     1097          double a00=aUnder[i+k*BLOCK]*multiplier;
     1098          double a01=aUnder[i+1+k*BLOCK]*multiplier;
     1099          t00 -= a00 * above[j + 0 + k * BLOCK];
     1100          t10 -= a00 * above[j + 1 + k * BLOCK];
     1101          t20 -= a00 * above[j + 2 + k * BLOCK];
     1102          t30 -= a00 * above[j + 3 + k * BLOCK];
     1103          t01 -= a01 * above[j + 0 + k * BLOCK];
     1104          t11 -= a01 * above[j + 1 + k * BLOCK];
     1105          t21 -= a01 * above[j + 2 + k * BLOCK];
     1106          t31 -= a01 * above[j + 3 + k * BLOCK];
     1107        }
     1108        aa[i+0*BLOCK]=t00;
     1109        aa[i+1*BLOCK]=t10;
     1110        aa[i+2*BLOCK]=t20;
     1111        aa[i+3*BLOCK]=t30;
     1112        aa[i+1+0*BLOCK]=t01;
     1113        aa[i+1+1*BLOCK]=t11;
     1114        aa[i+1+2*BLOCK]=t21;
     1115        aa[i+1+3*BLOCK]=t31;
     1116      }
     1117      if (odd) {
     1118        double t0=aa[n+0*BLOCK];
     1119        double t1=aa[n+1*BLOCK];
     1120        double t2=aa[n+2*BLOCK];
     1121        double t3=aa[n+3*BLOCK];
     1122        double a0;
     1123        for (k=0;k<BLOCK;k++) {
     1124          a0=aUnder[n+k*BLOCK]*work[k];
     1125          t0 -= a0 * above[j + 0 + k * BLOCK];
     1126          t1 -= a0 * above[j + 1 + k * BLOCK];
     1127          t2 -= a0 * above[j + 2 + k * BLOCK];
     1128          t3 -= a0 * above[j + 3 + k * BLOCK];
     1129        }
     1130        aa[n+0*BLOCK]=t0;
     1131        aa[n+1*BLOCK]=t1;
     1132        aa[n+2*BLOCK]=t2;
     1133        aa[n+3*BLOCK]=t3;
     1134      }
     1135    }
     1136  }
     1137#else
     1138  aa = aOther-BLOCK;
     1139  for (j = 0; j < BLOCK; j ++) {
     1140    aa+=BLOCK;
     1141    for (i=0;i< nUnder;i++) {
     1142      double t00=aa[i+0*BLOCK];
     1143      for (k=0;k<BLOCK;k++) {
     1144        double a00=aUnder[i+k*BLOCK]*work[k];
     1145        t00 -= a00 * above[j + k * BLOCK];
     1146      }
     1147      aa[i]=t00;
     1148    }
     1149  }
     1150#endif
    2091151}
    2101152/* Uses factorization to solve. */
     
    2121154ClpCholeskyDense::solve (double * region)
    2131155{
    214   int iRow,iColumn;
    215   for (iColumn=0;iColumn<numberRows_;iColumn++) {
    216     if (!rowsDropped_[iColumn]) {
    217       double value = region[iColumn];
    218       for (iRow=0;iRow<iColumn;iRow++)
    219         value -= region[iRow]*work_[iColumn+iRow*numberRows_];
    220       for (iRow=0;iRow<iColumn;iRow++)
    221         if (rowsDropped_[iRow])
    222           assert(!work_[iColumn+iRow*numberRows_]||!region[iRow]);
    223       region[iColumn]=value/work_[iColumn+iColumn*numberRows_];
     1156#ifdef CHOL_COMPARE
     1157  double * region2=NULL;
     1158  if (numberRows_<200) {
     1159    double * xx = sparseFactor_;
     1160    double * yy = diagonal_;
     1161    diagonal_ = sparseFactor_ + 40000;
     1162    sparseFactor_=diagonal_ + numberRows_;
     1163    region2 = ClpCopyOfArray(region,numberRows_);
     1164    int iRow,iColumn;
     1165    int addOffset=numberRows_-1;
     1166    double * work = sparseFactor_-1;
     1167    for (iColumn=0;iColumn<numberRows_;iColumn++) {
     1168      double value = region2[iColumn];
     1169      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
     1170        region2[iRow] -= value*work[iRow];
     1171      addOffset--;
     1172      work += addOffset;
     1173    }
     1174    for (iColumn=numberRows_-1;iColumn>=0;iColumn--) {
     1175      double value = region2[iColumn]*diagonal_[iColumn];
     1176      work -= addOffset;
     1177      addOffset++;
     1178      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
     1179        value -= region2[iRow]*work[iRow];
     1180      region2[iColumn]=value;
     1181    }
     1182    sparseFactor_=xx;
     1183    diagonal_=yy;
     1184  }
     1185#endif
     1186  //double * xx = sparseFactor_;
     1187  //double * yy = diagonal_;
     1188  //diagonal_ = sparseFactor_ + 40000;
     1189  //sparseFactor_=diagonal_ + numberRows_;
     1190  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
     1191  // later align on boundary
     1192  double * a = sparseFactor_+BLOCKSQ*numberBlocks;
     1193  int iBlock;
     1194  double * aa = a;
     1195  int iColumn;
     1196  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     1197    int nChunk;
     1198    int jBlock;
     1199    int iDo = iBlock*BLOCK;
     1200    int base=iDo;
     1201    if (iDo+BLOCK>numberRows_) {
     1202      nChunk=numberRows_-iDo;
    2241203    } else {
    225       region[iColumn]=0.0;
    226     }
    227   }
    228   double * work = work_ + numberRows_*numberRows_;
    229   for (iColumn=numberRows_-1;iColumn>=0;iColumn--) {
    230     work -= numberRows_;
    231     if (!rowsDropped_[iColumn]) {
    232       double value = region[iColumn];
    233       for (iRow=iColumn+1;iRow<numberRows_;iRow++)
    234         value -= region[iRow]*work[iRow];
    235       for (iRow=iColumn+1;iRow<numberRows_;iRow++)
    236         if (rowsDropped_[iRow])
    237           assert(!work[iRow]||!region[iRow]);
    238       region[iColumn]=value/work[iColumn];
     1204      nChunk=BLOCK;
     1205    }
     1206    solveF1(aa,nChunk,region+iDo);
     1207    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
     1208      base+=BLOCK;
     1209      aa+=BLOCKSQ;
     1210      if (base+BLOCK>numberRows_) {
     1211        nChunk=numberRows_-base;
     1212      } else {
     1213        nChunk=BLOCK;
     1214      }
     1215      solveF2(aa,nChunk,region+iDo,region+base);
     1216    }
     1217    aa+=BLOCKSQ;
     1218  }
     1219  // do diagonal outside
     1220  for (iColumn=0;iColumn<numberRows_;iColumn++)
     1221    region[iColumn] *= diagonal_[iColumn];
     1222  int offset=((numberBlocks*(numberBlocks+1))>>1);
     1223  aa = a+number_entries(offset-1);
     1224  int lBase=(numberBlocks-1)*BLOCK;
     1225  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
     1226    int nChunk;
     1227    int jBlock;
     1228    int triBase=iBlock*BLOCK;
     1229    int iBase=lBase;
     1230    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
     1231      if (iBase+BLOCK>numberRows_) {
     1232        nChunk=numberRows_-iBase;
     1233      } else {
     1234        nChunk=BLOCK;
     1235      }
     1236      solveB2(aa,nChunk,region+triBase,region+iBase);
     1237      iBase-=BLOCK;
     1238      aa-=BLOCKSQ;
     1239    }
     1240    if (triBase+BLOCK>numberRows_) {
     1241      nChunk=numberRows_-triBase;
    2391242    } else {
    240       region[iColumn]=0.0;
    241     }
    242   }
    243 }
     1243      nChunk=BLOCK;
     1244    }
     1245    solveB1(aa,nChunk,region+triBase);
     1246    aa-=BLOCKSQ;
     1247  }
     1248#ifdef CHOL_COMPARE
     1249  if (numberRows_<200) {
     1250    for (int i=0;i<numberRows_;i++) {
     1251      assert(fabs(region[i]-region2[i])<1.0e-3);
     1252    }
     1253    delete [] region2;
     1254  }
     1255#endif
     1256}
     1257// Forward part of solve 1
     1258void
     1259ClpCholeskyDense::solveF1(double * a,int n,double * region)
     1260{
     1261  int j, k;
     1262  double t00;
     1263  for (j = 0; j < n; j ++) {
     1264    t00 = region[j];
     1265    for (k = 0; k < j; ++k) {
     1266      t00 -= region[k] * a[j + k * BLOCK];
     1267    }
     1268    //t00*=a[j + j * BLOCK];
     1269    region[j] = t00;
     1270  }
     1271}
     1272// Forward part of solve 2
     1273void
     1274ClpCholeskyDense::solveF2(double * a,int n,double * region, double * region2)
     1275{
     1276  int j, k;
     1277#if BLOCK==16
     1278  if (n==BLOCK) {
     1279    for (k = 0; k < BLOCK; k+=4) {
     1280      double t0 = region2[0];
     1281      double t1 = region2[1];
     1282      double t2 = region2[2];
     1283      double t3 = region2[3];
     1284      t0 -= region[0] * a[0 + 0 * BLOCK];
     1285      t1 -= region[0] * a[1 + 0 * BLOCK];
     1286      t2 -= region[0] * a[2 + 0 * BLOCK];
     1287      t3 -= region[0] * a[3 + 0 * BLOCK];
     1288
     1289      t0 -= region[1] * a[0 + 1 * BLOCK];
     1290      t1 -= region[1] * a[1 + 1 * BLOCK];
     1291      t2 -= region[1] * a[2 + 1 * BLOCK];
     1292      t3 -= region[1] * a[3 + 1 * BLOCK];
     1293
     1294      t0 -= region[2] * a[0 + 2 * BLOCK];
     1295      t1 -= region[2] * a[1 + 2 * BLOCK];
     1296      t2 -= region[2] * a[2 + 2 * BLOCK];
     1297      t3 -= region[2] * a[3 + 2 * BLOCK];
     1298
     1299      t0 -= region[3] * a[0 + 3 * BLOCK];
     1300      t1 -= region[3] * a[1 + 3 * BLOCK];
     1301      t2 -= region[3] * a[2 + 3 * BLOCK];
     1302      t3 -= region[3] * a[3 + 3 * BLOCK];
     1303
     1304      t0 -= region[4] * a[0 + 4 * BLOCK];
     1305      t1 -= region[4] * a[1 + 4 * BLOCK];
     1306      t2 -= region[4] * a[2 + 4 * BLOCK];
     1307      t3 -= region[4] * a[3 + 4 * BLOCK];
     1308
     1309      t0 -= region[5] * a[0 + 5 * BLOCK];
     1310      t1 -= region[5] * a[1 + 5 * BLOCK];
     1311      t2 -= region[5] * a[2 + 5 * BLOCK];
     1312      t3 -= region[5] * a[3 + 5 * BLOCK];
     1313
     1314      t0 -= region[6] * a[0 + 6 * BLOCK];
     1315      t1 -= region[6] * a[1 + 6 * BLOCK];
     1316      t2 -= region[6] * a[2 + 6 * BLOCK];
     1317      t3 -= region[6] * a[3 + 6 * BLOCK];
     1318
     1319      t0 -= region[7] * a[0 + 7 * BLOCK];
     1320      t1 -= region[7] * a[1 + 7 * BLOCK];
     1321      t2 -= region[7] * a[2 + 7 * BLOCK];
     1322      t3 -= region[7] * a[3 + 7 * BLOCK];
     1323#if BLOCK>8
     1324      t0 -= region[8] * a[0 + 8 * BLOCK];
     1325      t1 -= region[8] * a[1 + 8 * BLOCK];
     1326      t2 -= region[8] * a[2 + 8 * BLOCK];
     1327      t3 -= region[8] * a[3 + 8 * BLOCK];
     1328
     1329      t0 -= region[9] * a[0 + 9 * BLOCK];
     1330      t1 -= region[9] * a[1 + 9 * BLOCK];
     1331      t2 -= region[9] * a[2 + 9 * BLOCK];
     1332      t3 -= region[9] * a[3 + 9 * BLOCK];
     1333
     1334      t0 -= region[10] * a[0 + 10 * BLOCK];
     1335      t1 -= region[10] * a[1 + 10 * BLOCK];
     1336      t2 -= region[10] * a[2 + 10 * BLOCK];
     1337      t3 -= region[10] * a[3 + 10 * BLOCK];
     1338
     1339      t0 -= region[11] * a[0 + 11 * BLOCK];
     1340      t1 -= region[11] * a[1 + 11 * BLOCK];
     1341      t2 -= region[11] * a[2 + 11 * BLOCK];
     1342      t3 -= region[11] * a[3 + 11 * BLOCK];
     1343
     1344      t0 -= region[12] * a[0 + 12 * BLOCK];
     1345      t1 -= region[12] * a[1 + 12 * BLOCK];
     1346      t2 -= region[12] * a[2 + 12 * BLOCK];
     1347      t3 -= region[12] * a[3 + 12 * BLOCK];
     1348
     1349      t0 -= region[13] * a[0 + 13 * BLOCK];
     1350      t1 -= region[13] * a[1 + 13 * BLOCK];
     1351      t2 -= region[13] * a[2 + 13 * BLOCK];
     1352      t3 -= region[13] * a[3 + 13 * BLOCK];
     1353
     1354      t0 -= region[14] * a[0 + 14 * BLOCK];
     1355      t1 -= region[14] * a[1 + 14 * BLOCK];
     1356      t2 -= region[14] * a[2 + 14 * BLOCK];
     1357      t3 -= region[14] * a[3 + 14 * BLOCK];
     1358
     1359      t0 -= region[15] * a[0 + 15 * BLOCK];
     1360      t1 -= region[15] * a[1 + 15 * BLOCK];
     1361      t2 -= region[15] * a[2 + 15 * BLOCK];
     1362      t3 -= region[15] * a[3 + 15 * BLOCK];
     1363#endif
     1364      region2[0] = t0;
     1365      region2[1] = t1;
     1366      region2[2] = t2;
     1367      region2[3] = t3;
     1368      region2+=4;
     1369      a+=4;
     1370    }
     1371  } else {
     1372#endif
     1373    for (k = 0; k < n; ++k) {
     1374      double t00 = region2[k];
     1375      for (j = 0; j < BLOCK; j ++) {
     1376        t00 -= region[j] * a[k + j * BLOCK];
     1377      }
     1378      region2[k] = t00;
     1379    }
     1380#if BLOCK==16
     1381  }
     1382#endif
     1383}
     1384// Backward part of solve 1
     1385void
     1386ClpCholeskyDense::solveB1(double * a,int n,double * region)
     1387{
     1388  int j, k;
     1389  double t00;
     1390  for (j = n-1; j >=0; j --) {
     1391    t00 = region[j];
     1392    for (k = j+1; k < n; ++k) {
     1393      t00 -= region[k] * a[k + j * BLOCK];
     1394    }
     1395    //t00*=a[j + j * BLOCK];
     1396    region[j] = t00;
     1397  }
     1398}
     1399// Backward part of solve 2
     1400void
     1401ClpCholeskyDense::solveB2(double * a,int n,double * region, double * region2)
     1402{
     1403  int j, k;
     1404#if BLOCK==16
     1405  if (n==BLOCK) {
     1406    for (j = 0; j < BLOCK; j +=4) {
     1407      double t0 = region[0];
     1408      double t1 = region[1];
     1409      double t2 = region[2];
     1410      double t3 = region[3];
     1411      t0 -= region2[0] * a[0 + 0*BLOCK];
     1412      t1 -= region2[0] * a[0 + 1*BLOCK];
     1413      t2 -= region2[0] * a[0 + 2*BLOCK];
     1414      t3 -= region2[0] * a[0 + 3*BLOCK];
     1415
     1416      t0 -= region2[1] * a[1 + 0*BLOCK];
     1417      t1 -= region2[1] * a[1 + 1*BLOCK];
     1418      t2 -= region2[1] * a[1 + 2*BLOCK];
     1419      t3 -= region2[1] * a[1 + 3*BLOCK];
     1420
     1421      t0 -= region2[2] * a[2 + 0*BLOCK];
     1422      t1 -= region2[2] * a[2 + 1*BLOCK];
     1423      t2 -= region2[2] * a[2 + 2*BLOCK];
     1424      t3 -= region2[2] * a[2 + 3*BLOCK];
     1425
     1426      t0 -= region2[3] * a[3 + 0*BLOCK];
     1427      t1 -= region2[3] * a[3 + 1*BLOCK];
     1428      t2 -= region2[3] * a[3 + 2*BLOCK];
     1429      t3 -= region2[3] * a[3 + 3*BLOCK];
     1430
     1431      t0 -= region2[4] * a[4 + 0*BLOCK];
     1432      t1 -= region2[4] * a[4 + 1*BLOCK];
     1433      t2 -= region2[4] * a[4 + 2*BLOCK];
     1434      t3 -= region2[4] * a[4 + 3*BLOCK];
     1435
     1436      t0 -= region2[5] * a[5 + 0*BLOCK];
     1437      t1 -= region2[5] * a[5 + 1*BLOCK];
     1438      t2 -= region2[5] * a[5 + 2*BLOCK];
     1439      t3 -= region2[5] * a[5 + 3*BLOCK];
     1440
     1441      t0 -= region2[6] * a[6 + 0*BLOCK];
     1442      t1 -= region2[6] * a[6 + 1*BLOCK];
     1443      t2 -= region2[6] * a[6 + 2*BLOCK];
     1444      t3 -= region2[6] * a[6 + 3*BLOCK];
     1445
     1446      t0 -= region2[7] * a[7 + 0*BLOCK];
     1447      t1 -= region2[7] * a[7 + 1*BLOCK];
     1448      t2 -= region2[7] * a[7 + 2*BLOCK];
     1449      t3 -= region2[7] * a[7 + 3*BLOCK];
     1450#if BLOCK>8
     1451
     1452      t0 -= region2[8] * a[8 + 0*BLOCK];
     1453      t1 -= region2[8] * a[8 + 1*BLOCK];
     1454      t2 -= region2[8] * a[8 + 2*BLOCK];
     1455      t3 -= region2[8] * a[8 + 3*BLOCK];
     1456
     1457      t0 -= region2[9] * a[9 + 0*BLOCK];
     1458      t1 -= region2[9] * a[9 + 1*BLOCK];
     1459      t2 -= region2[9] * a[9 + 2*BLOCK];
     1460      t3 -= region2[9] * a[9 + 3*BLOCK];
     1461
     1462      t0 -= region2[10] * a[10 + 0*BLOCK];
     1463      t1 -= region2[10] * a[10 + 1*BLOCK];
     1464      t2 -= region2[10] * a[10 + 2*BLOCK];
     1465      t3 -= region2[10] * a[10 + 3*BLOCK];
     1466
     1467      t0 -= region2[11] * a[11 + 0*BLOCK];
     1468      t1 -= region2[11] * a[11 + 1*BLOCK];
     1469      t2 -= region2[11] * a[11 + 2*BLOCK];
     1470      t3 -= region2[11] * a[11 + 3*BLOCK];
     1471
     1472      t0 -= region2[12] * a[12 + 0*BLOCK];
     1473      t1 -= region2[12] * a[12 + 1*BLOCK];
     1474      t2 -= region2[12] * a[12 + 2*BLOCK];
     1475      t3 -= region2[12] * a[12 + 3*BLOCK];
     1476
     1477      t0 -= region2[13] * a[13 + 0*BLOCK];
     1478      t1 -= region2[13] * a[13 + 1*BLOCK];
     1479      t2 -= region2[13] * a[13 + 2*BLOCK];
     1480      t3 -= region2[13] * a[13 + 3*BLOCK];
     1481
     1482      t0 -= region2[14] * a[14 + 0*BLOCK];
     1483      t1 -= region2[14] * a[14 + 1*BLOCK];
     1484      t2 -= region2[14] * a[14 + 2*BLOCK];
     1485      t3 -= region2[14] * a[14 + 3*BLOCK];
     1486
     1487      t0 -= region2[15] * a[15 + 0*BLOCK];
     1488      t1 -= region2[15] * a[15 + 1*BLOCK];
     1489      t2 -= region2[15] * a[15 + 2*BLOCK];
     1490      t3 -= region2[15] * a[15 + 3*BLOCK];
     1491#endif
     1492      region[0] = t0;
     1493      region[1] = t1;
     1494      region[2] = t2;
     1495      region[3] = t3;
     1496      a+=4*BLOCK;
     1497      region+=4;
     1498    }
     1499  } else {
     1500#endif
     1501    for (j = 0; j < BLOCK; j ++) {
     1502      double t00 = region[j];
     1503      for (k = 0; k < n; ++k) {
     1504        t00 -= region2[k] * a[k + j * BLOCK];
     1505      }
     1506      region[j] = t00;
     1507    }
     1508#if BLOCK==16
     1509  }
     1510#endif
     1511}
  • trunk/ClpCholeskyTaucs.cpp

    r399 r414  
    2424    matrix_(NULL),
    2525    factorization_(NULL),
    26     sparseFactor_(NULL),
    27     choleskyStart_(NULL),
    28     choleskyRow_(NULL),
    29     sizeFactor_(0),
    30     rowCopy_(NULL)
     26    sparseFactorT_(NULL),
     27    choleskyStartT_(NULL),
     28    choleskyRowT_(NULL),
     29    sizeFactorT_(0),
     30    rowCopyT_(NULL)
    3131{
    3232  type_=13;
     
    4242  // For Taucs stuff is done by malloc
    4343  matrix_ = rhs.matrix_;
    44   sizeFactor_=rhs.sizeFactor_;
     44  sizeFactorT_=rhs.sizeFactorT_;
    4545  if (matrix_) {
    46     choleskyStart_ = (int *) malloc((numberRows_+1)*sizeof(int));
    47     memcpy(choleskyStart_,rhs.choleskyStart_,(numberRows_+1)*sizeof(int));
    48     choleskyRow_ = (int *) malloc(sizeFactor_*sizeof(int));
    49     memcpy(choleskyRow_,rhs.choleskyRow_,sizeFactor_*sizeof(int));
    50     sparseFactor_ = (double *) malloc(sizeFactor_*sizeof(double));
    51     memcpy(sparseFactor_,rhs.sparseFactor_,sizeFactor_*sizeof(double));
    52     matrix_->colptr = choleskyStart_;
    53     matrix_->rowind = choleskyRow_;
    54     matrix_->values.d = sparseFactor_;
     46    choleskyStartT_ = (int *) malloc((numberRows_+1)*sizeof(int));
     47    memcpy(choleskyStartT_,rhs.choleskyStartT_,(numberRows_+1)*sizeof(int));
     48    choleskyRowT_ = (int *) malloc(sizeFactorT_*sizeof(int));
     49    memcpy(choleskyRowT_,rhs.choleskyRowT_,sizeFactorT_*sizeof(int));
     50    sparseFactorT_ = (double *) malloc(sizeFactorT_*sizeof(double));
     51    memcpy(sparseFactorT_,rhs.sparseFactorT_,sizeFactorT_*sizeof(double));
     52    matrix_->colptr = choleskyStartT_;
     53    matrix_->rowind = choleskyRowT_;
     54    matrix_->values.d = sparseFactorT_;
    5555  } else {
    56     sparseFactor_=NULL;
    57     choleskyStart_=NULL;
    58     choleskyRow_=NULL;
     56    sparseFactorT_=NULL;
     57    choleskyStartT_=NULL;
     58    choleskyRowT_=NULL;
    5959  }
    6060  factorization_=NULL,
    61   rowCopy_ = rhs.rowCopy_->clone();
     61  rowCopyT_ = rhs.rowCopyT_->clone();
    6262}
    6363
     
    7171  if (factorization_)
    7272    taucs_supernodal_factor_free(factorization_);
    73   delete rowCopy_;
     73  delete rowCopyT_;
    7474}
    7575
     
    8686      taucs_supernodal_factor_free(factorization_);
    8787    factorization_=NULL;
    88     sizeFactor_=rhs.sizeFactor_;
     88    sizeFactorT_=rhs.sizeFactorT_;
    8989    matrix_ = rhs.matrix_;
    9090    if (matrix_) {
    91       choleskyStart_ = (int *) malloc((numberRows_+1)*sizeof(int));
    92       memcpy(choleskyStart_,rhs.choleskyStart_,(numberRows_+1)*sizeof(int));
    93       choleskyRow_ = (int *) malloc(sizeFactor_*sizeof(int));
    94       memcpy(choleskyRow_,rhs.choleskyRow_,sizeFactor_*sizeof(int));
    95       sparseFactor_ = (double *) malloc(sizeFactor_*sizeof(double));
    96       memcpy(sparseFactor_,rhs.sparseFactor_,sizeFactor_*sizeof(double));
    97       matrix_->colptr = choleskyStart_;
    98       matrix_->rowind = choleskyRow_;
    99       matrix_->values.d = sparseFactor_;
     91      choleskyStartT_ = (int *) malloc((numberRows_+1)*sizeof(int));
     92      memcpy(choleskyStartT_,rhs.choleskyStartT_,(numberRows_+1)*sizeof(int));
     93      choleskyRowT_ = (int *) malloc(sizeFactorT_*sizeof(int));
     94      memcpy(choleskyRowT_,rhs.choleskyRowT_,sizeFactorT_*sizeof(int));
     95      sparseFactorT_ = (double *) malloc(sizeFactorT_*sizeof(double));
     96      memcpy(sparseFactorT_,rhs.sparseFactorT_,sizeFactorT_*sizeof(double));
     97      matrix_->colptr = choleskyStartT_;
     98      matrix_->rowind = choleskyRowT_;
     99      matrix_->values.d = sparseFactorT_;
    100100    } else {
    101       sparseFactor_=NULL;
    102       choleskyStart_=NULL;
    103       choleskyRow_=NULL;
     101      sparseFactorT_=NULL;
     102      choleskyStartT_=NULL;
     103      choleskyRowT_=NULL;
    104104    }
    105     delete rowCopy_;
    106     rowCopy_ = rhs.rowCopy_->clone();
     105    delete rowCopyT_;
     106    rowCopyT_ = rhs.rowCopyT_->clone();
    107107  }
    108108  return *this;
     
    124124  numberRowsDropped_=0;
    125125  model_=model;
    126   rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     126  rowCopyT_ = model->clpMatrix()->reverseOrderedCopy();
    127127  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    128128  const int * columnLength = model_->clpMatrix()->getVectorLengths();
    129129  const int * row = model_->clpMatrix()->getIndices();
    130   const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    131   const int * rowLength = rowCopy_->getVectorLengths();
    132   const int * column = rowCopy_->getIndices();
     130  const CoinBigIndex * rowStart = rowCopyT_->getVectorStarts();
     131  const int * rowLength = rowCopyT_->getVectorLengths();
     132  const int * column = rowCopyT_->getIndices();
    133133  // We need two arrays for counts
    134134  int * which = new int [numberRows_];
     
    136136  CoinZeroN(used,numberRows_);
    137137  int iRow;
    138   sizeFactor_=0;
     138  sizeFactorT_=0;
    139139  for (iRow=0;iRow<numberRows_;iRow++) {
    140140    int number=1;
     
    159159        }
    160160      }
    161       sizeFactor_ += number;
     161      sizeFactorT_ += number;
    162162      int j;
    163163      for (j=0;j<number;j++)
     
    167167  delete [] which;
    168168  // Now we have size - create arrays and fill in
    169   matrix_ = taucs_ccs_create(numberRows_,numberRows_,sizeFactor_,
     169  matrix_ = taucs_ccs_create(numberRows_,numberRows_,sizeFactorT_,
    170170                             TAUCS_DOUBLE|TAUCS_SYMMETRIC|TAUCS_LOWER);
    171171  if (!matrix_)
    172172    return 1;
    173173  // Space for starts
    174   choleskyStart_ = matrix_->colptr;
    175   choleskyRow_ = matrix_->rowind;
    176   sparseFactor_ = matrix_->values.d;
    177   sizeFactor_=0;
    178   which = choleskyRow_;
     174  choleskyStartT_ = matrix_->colptr;
     175  choleskyRowT_ = matrix_->rowind;
     176  sparseFactorT_ = matrix_->values.d;
     177  sizeFactorT_=0;
     178  which = choleskyRowT_;
    179179  for (iRow=0;iRow<numberRows_;iRow++) {
    180180    int number=1;
     
    182182    which[0]=iRow;
    183183    used[iRow]=1;
    184     choleskyStart_[iRow]=sizeFactor_;
     184    choleskyStartT_[iRow]=sizeFactorT_;
    185185    if (!rowsDropped_[iRow]) {
    186186      CoinBigIndex startRow=rowStart[iRow];
     
    200200        }
    201201      }
    202       sizeFactor_ += number;
     202      sizeFactorT_ += number;
    203203      int j;
    204204      for (j=0;j<number;j++)
     
    210210    }
    211211  }
    212   choleskyStart_[numberRows_]=sizeFactor_;
     212  choleskyStartT_[numberRows_]=sizeFactorT_;
    213213  delete [] used;
    214   permuteIn_ = new int [numberRows_];
    215   permuteOut_ = new int[numberRows_];
     214  permuteInverse_ = new int [numberRows_];
     215  permute_ = new int[numberRows_];
    216216  int * perm, *invp;
    217217  // There seem to be bugs in ordering if model too small
     
    220220  else
    221221    taucs_ccs_order(matrix_,&perm,&invp,(const char *) "identity");
    222   memcpy(permuteIn_,perm,numberRows_*sizeof(int));
     222  memcpy(permuteInverse_,perm,numberRows_*sizeof(int));
    223223  free(perm);
    224   memcpy(permuteOut_,invp,numberRows_*sizeof(int));
     224  memcpy(permute_,invp,numberRows_*sizeof(int));
    225225  free(invp);
    226226  // need to permute
    227   taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_,permuteIn_,permuteOut_);
     227  taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_,permuteInverse_,permute_);
    228228  // symbolic
    229229  factorization_ = taucs_ccs_factor_llt_symbolic(permuted);
     
    239239  const int * row = model_->clpMatrix()->getIndices();
    240240  const double * element = model_->clpMatrix()->getElements();
    241   const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
    242   const int * rowLength = rowCopy_->getVectorLengths();
    243   const int * column = rowCopy_->getIndices();
    244   const double * elementByRow = rowCopy_->getElements();
     241  const CoinBigIndex * rowStart = rowCopyT_->getVectorStarts();
     242  const int * rowLength = rowCopyT_->getVectorLengths();
     243  const int * column = rowCopyT_->getIndices();
     244  const double * elementByRow = rowCopyT_->getElements();
    245245  int numberColumns=model_->clpMatrix()->getNumCols();
    246246  int iRow;
     
    260260  }
    261261  for (iRow=0;iRow<numberRows_;iRow++) {
    262     double * put = sparseFactor_+choleskyStart_[iRow];
    263     int * which = choleskyRow_+choleskyStart_[iRow];
    264     int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
     262    double * put = sparseFactorT_+choleskyStartT_[iRow];
     263    int * which = choleskyRowT_+choleskyStartT_[iRow];
     264    int number = choleskyStartT_[iRow+1]-choleskyStartT_[iRow];
    265265    if (!rowLength[iRow])
    266266      rowsDropped_[iRow]=1;
     
    298298  }
    299299  //check sizes
    300   double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
     300  double largest2=maximumAbsElement(sparseFactorT_,sizeFactorT_);
    301301  largest2*=1.0e-19;
    302302  largest = CoinMin(largest2,1.0e-11);
     
    307307    rowsDropped[iRow]=dropped;
    308308    if (!dropped) {
    309       CoinBigIndex start = choleskyStart_[iRow];
    310       double diagonal = sparseFactor_[start];
     309      CoinBigIndex start = choleskyStartT_[iRow];
     310      double diagonal = sparseFactorT_[start];
    311311      if (diagonal>largest2) {
    312         sparseFactor_[start]=diagonal+perturbation;
     312        sparseFactorT_[start]=diagonal+perturbation;
    313313      } else {
    314         sparseFactor_[start]=diagonal+perturbation;
     314        sparseFactorT_[start]=diagonal+perturbation;
    315315        rowsDropped[iRow]=2;
    316316        numberDroppedBefore++;
     
    320320  taucs_supernodal_factor_free_numeric(factorization_);
    321321  // need to permute
    322   taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_,permuteIn_,permuteOut_);
     322  taucs_ccs_matrix * permuted = taucs_ccs_permute_symmetrically(matrix_,permuteInverse_,permute_);
    323323  int rCode=taucs_ccs_factor_llt_numeric(permuted,factorization_);
    324324  taucs_ccs_free(permuted);
     
    391391  double * in = new double[numberRows_];
    392392  double * out = new double[numberRows_];
    393   taucs_vec_permute(numberRows_,TAUCS_DOUBLE,region,in,permuteIn_);
     393  taucs_vec_permute(numberRows_,TAUCS_DOUBLE,region,in,permuteInverse_);
    394394  int rCode=taucs_supernodal_solve_llt(factorization_,out,in);
    395395  if (rCode)
    396396    printf("return code of %d from solve\n",rCode);
    397   taucs_vec_permute(numberRows_,TAUCS_DOUBLE,out,region,permuteOut_);
     397  taucs_vec_permute(numberRows_,TAUCS_DOUBLE,out,region,permute_);
    398398  delete [] out;
    399399  delete [] in;
  • trunk/ClpCholeskyWssmp.cpp

    r399 r414  
    2222//-------------------------------------------------------------------
    2323ClpCholeskyWssmp::ClpCholeskyWssmp (int denseThreshold)
    24   : ClpCholeskyBase(),
    25     sparseFactor_(NULL),
    26     choleskyStart_(NULL),
    27     choleskyRow_(NULL),
    28     sizeFactor_(0),
    29     rowCopy_(NULL),
    30     whichDense_(NULL),
    31     denseColumn_(NULL),
    32     dense_(NULL),
    33     denseThreshold_(denseThreshold)
     24  : ClpCholeskyBase(denseThreshold)
    3425{
    3526  type_=12;
    36   memset(integerParameters_,0,64*sizeof(int));
    37   memset(doubleParameters_,0,64*sizeof(double));
    3827}
    3928
     
    4433: ClpCholeskyBase(rhs)
    4534{
    46   type_=rhs.type_;
    47   sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
    48   choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
    49   choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
    50   sizeFactor_=rhs.sizeFactor_;
    51   memcpy(integerParameters_,rhs.integerParameters_,64*sizeof(int));
    52   memcpy(doubleParameters_,rhs.doubleParameters_,64*sizeof(double));
    53   rowCopy_ = rhs.rowCopy_->clone();
    54   whichDense_ = NULL;
    55   denseColumn_=NULL;
    56   dense_=NULL;
    57   denseThreshold_ = rhs.denseThreshold_;
    5835}
    5936
     
    6441ClpCholeskyWssmp::~ClpCholeskyWssmp ()
    6542{
    66   delete [] sparseFactor_;
    67   delete [] choleskyStart_;
    68   delete [] choleskyRow_;
    69   delete rowCopy_;
    70   delete [] whichDense_;
    71   delete [] denseColumn_;
    72   delete dense_;
    7343}
    7444
     
    8151  if (this != &rhs) {
    8252    ClpCholeskyBase::operator=(rhs);
    83     delete [] sparseFactor_;
    84     delete [] choleskyStart_;
    85     delete [] choleskyRow_;
    86     delete [] whichDense_;
    87     delete [] denseColumn_;
    88     delete dense_;
    89     sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
    90     choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
    91     choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
    92     sizeFactor_=rhs.sizeFactor_;
    93     delete rowCopy_;
    94     rowCopy_ = rhs.rowCopy_->clone();
    95     whichDense_ = NULL;
    96     denseColumn_=NULL;
    97     dense_=NULL;
    98     denseThreshold_ = rhs.denseThreshold_;
    9953  }
    10054  return *this;
     
    222176      // dense cholesky
    223177      dense_ = new ClpCholeskyDense();
    224       dense_->reserveSpace(numberDense);
     178      dense_->reserveSpace(NULL,numberDense);
    225179      printf("Taking %d columns as dense\n",numberDense);
    226180    }
     
    318272  choleskyStart_[numberRows_]=sizeFactor_;
    319273  delete [] used;
    320   permuteIn_ = new int [numberRows_];
    321   permuteOut_ = new int[numberRows_];
     274  permuteInverse_ = new int [numberRows_];
     275  permute_ = new int[numberRows_];
    322276  integerParameters_[0]=0;
    323277  int i0=0;
     
    327281#endif
    328282  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    329          NULL,permuteOut_,permuteIn_,0,&numberRows_,&i1,
     283         NULL,permute_,permuteInverse_,0,&numberRows_,&i1,
    330284         NULL,&i0,NULL,integerParameters_,doubleParameters_);
    331285  integerParameters_[1]=1;//order and symbolic
     
    342296  doubleParameters_[11]=1.0e-15;
    343297  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    344          NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
     298         NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
    345299         NULL,&i0,NULL,integerParameters_,doubleParameters_);
    346300  //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
     
    352306  if (!integerParameters_[23]) {
    353307    for (int iRow=0;iRow<numberRows_;iRow++) {
    354       permuteIn_[iRow]=iRow;
    355       permuteOut_[iRow]=iRow;
     308      permuteInverse_[iRow]=iRow;
     309      permute_[iRow]=iRow;
    356310    }
    357311    std::cout<<"wssmp says no elements - fully dense? - switching to dense"<<std::endl;
     
    360314    integerParameters_[7]=1; // no permute
    361315    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    362           NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
     316          NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
    363317          NULL,&i0,NULL,integerParameters_,doubleParameters_);
    364318    std::cout<<integerParameters_[23]<<" elements in dense Cholesky"<<std::endl;
    365319  }
     320  return 0;
     321}
     322/* Does Symbolic factorization given permutation.
     323   This is called immediately after order.  If user provides this then
     324   user must provide factorize and solve.  Otherwise the default factorization is used
     325   returns non-zero if not enough memory */
     326int
     327ClpCholeskyWssmp::symbolic()
     328{
    366329  return 0;
    367330}
     
    399362  }
    400363  if (whichDense_) {
    401     longDouble * denseBlob = dense_->aMatrix();
    402     CoinZeroN(denseBlob,numberDense*numberDense);
     364    double * denseDiagonal = dense_->diagonal();
    403365    double * dense = denseColumn_;
    404     longDouble * blob = denseBlob;
     366    int iDense=0;
    405367    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    406368      if (whichDense_[iColumn]) {
    407         blob[0]=1.0/diagonal[iColumn];
     369        denseDiagonal[iDense++]=1.0/diagonal[iColumn];
    408370        CoinZeroN(dense,numberRows_);
    409371        CoinBigIndex start=columnStart[iColumn];
     
    414376        }
    415377        dense += numberRows_;
    416         blob += numberDense+1;
    417378      }
    418379    }
     
    494455    integerParameters_[9]=0;
    495456  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    496         NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
     457        NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
    497458        NULL,&i0,rowsDropped,integerParameters_,doubleParameters_);
    498459  //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
     
    527488    integerParameters_[2]=4;
    528489    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    529        NULL,permuteOut_,permuteIn_,denseColumn_,&numberRows_,&numberDense,
     490       NULL,permute_,permuteInverse_,denseColumn_,&numberRows_,&numberDense,
    530491       NULL,&i0,NULL,integerParameters_,doubleParameters_);
    531492    integerParameters_[29]=0;
    532493    dense_->resetRowsDropped();
    533494    longDouble * denseBlob = dense_->aMatrix();
     495    double * denseDiagonal = dense_->diagonal();
    534496    // Update dense matrix
    535497    for (i=0;i<numberDense;i++) {
    536498      const double * a = denseColumn_+i*numberRows_;
    537       for (int j=0;j<numberDense;j++) {
    538         double value = denseBlob[i+j*numberDense];
     499      // do diagonal
     500      double value = denseDiagonal[i];
     501      const double * b = denseColumn_+i*numberRows_;
     502      for (int k=0;k<numberRows_;k++)
     503        value += a[k]*b[k];
     504      denseDiagonal[i]=value;
     505      for (int j=i+1;j<numberDense;j++) {
     506        double value = 0.0;
    539507        const double * b = denseColumn_+j*numberRows_;
    540508        for (int k=0;k<numberRows_;k++)
    541509          value += a[k]*b[k];
    542         denseBlob[i+j*numberDense]=value;
     510        *denseBlob=value;
     511        denseBlob++;
    543512      }
    544513    }
     
    616585  if (!whichDense_) {
    617586    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    618           NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
     587          NULL,permute_,permuteInverse_,region,&numberRows_,&i1,
    619588          NULL,&i0,NULL,integerParameters_,doubleParameters_);
    620589  } else {
     
    622591    integerParameters_[29]=1;
    623592    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    624           NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
     593          NULL,permute_,permuteInverse_,region,&numberRows_,&i1,
    625594          NULL,&i0,NULL,integerParameters_,doubleParameters_);
    626595    // do change;
     
    648617    integerParameters_[1]=4;
    649618    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    650           NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
     619          NULL,permute_,permuteInverse_,region,&numberRows_,&i1,
    651620          NULL,&i0,NULL,integerParameters_,doubleParameters_);
    652621    integerParameters_[29]=0;
  • trunk/ClpCholeskyWssmpKKT.cpp

    r399 r414  
    2222//-------------------------------------------------------------------
    2323ClpCholeskyWssmpKKT::ClpCholeskyWssmpKKT (int denseThreshold)
    24   : ClpCholeskyBase(),
    25     sparseFactor_(NULL),
    26     choleskyStart_(NULL),
    27     choleskyRow_(NULL),
    28     sizeFactor_(0),
    29     denseThreshold_(denseThreshold)
     24  : ClpCholeskyBase(denseThreshold)
    3025{
    3126  type_=21;
    32   memset(integerParameters_,0,64*sizeof(int));
    33   memset(doubleParameters_,0,64*sizeof(double));
    3427}
    3528
     
    4033: ClpCholeskyBase(rhs)
    4134{
    42   type_=rhs.type_;
    43   sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
    44   choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
    45   choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
    46   sizeFactor_=rhs.sizeFactor_;
    47   memcpy(integerParameters_,rhs.integerParameters_,64*sizeof(int));
    48   memcpy(doubleParameters_,rhs.doubleParameters_,64*sizeof(double));
    49   denseThreshold_ = rhs.denseThreshold_;
    5035}
    5136
     
    5641ClpCholeskyWssmpKKT::~ClpCholeskyWssmpKKT ()
    5742{
    58   delete [] sparseFactor_;
    59   delete [] choleskyStart_;
    60   delete [] choleskyRow_;
    6143}
    6244
     
    6951  if (this != &rhs) {
    7052    ClpCholeskyBase::operator=(rhs);
    71     delete [] sparseFactor_;
    72     delete [] choleskyStart_;
    73     delete [] choleskyRow_;
    74     sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
    75     choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
    76     choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
    77     sizeFactor_=rhs.sizeFactor_;
    78     denseThreshold_ = rhs.denseThreshold_;
    7953  }
    8054  return *this;
     
    233207  }
    234208  choleskyStart_[numberRows_]=sizeFactor_;
    235   permuteIn_ = new int [numberRows_];
    236   permuteOut_ = new int[numberRows_];
     209  permuteInverse_ = new int [numberRows_];
     210  permute_ = new int[numberRows_];
    237211  integerParameters_[0]=0;
    238212  int i0=0;
     
    242216#endif
    243217  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    244          NULL,permuteOut_,permuteIn_,0,&numberRows_,&i1,
     218         NULL,permute_,permuteInverse_,0,&numberRows_,&i1,
    245219         NULL,&i0,NULL,integerParameters_,doubleParameters_);
    246220  integerParameters_[1]=1;//order and symbolic
     
    259233  integerParameters_[1]=2;//just symbolic
    260234  for (int iRow=0;iRow<numberRows_;iRow++) {
    261     permuteIn_[iRow]=iRow;
    262     permuteOut_[iRow]=iRow;
     235    permuteInverse_[iRow]=iRow;
     236    permute_[iRow]=iRow;
    263237  }
    264238#endif
    265239  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    266          NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
     240         NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
    267241         NULL,&i0,NULL,integerParameters_,doubleParameters_);
    268242  //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
     
    274248  if (!integerParameters_[23]) {
    275249    for (int iRow=0;iRow<numberRows_;iRow++) {
    276       permuteIn_[iRow]=iRow;
    277       permuteOut_[iRow]=iRow;
     250      permuteInverse_[iRow]=iRow;
     251      permute_[iRow]=iRow;
    278252    }
    279253    std::cout<<"wssmp says no elements - fully dense? - switching to dense"<<std::endl;
     
    282256    integerParameters_[7]=1; // no permute
    283257    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    284           NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
     258          NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
    285259          NULL,&i0,NULL,integerParameters_,doubleParameters_);
    286260    std::cout<<integerParameters_[23]<<" elements in dense Cholesky"<<std::endl;
    287261  }
     262  return 0;
     263}
     264/* Does Symbolic factorization given permutation.
     265   This is called immediately after order.  If user provides this then
     266   user must provide factorize and solve.  Otherwise the default factorization is used
     267   returns non-zero if not enough memory */
     268int
     269ClpCholeskyWssmpKKT::symbolic()
     270{
    288271  return 0;
    289272}
     
    435418  CoinZeroN(rowsDropped2,numberRows_);
    436419  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    437         NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
     420        NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
    438421        NULL,&i0,rowsDropped2,integerParameters_,doubleParameters_);
    439422   //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
     
    513496#endif
    514497  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
    515         NULL,permuteOut_,permuteIn_,array,&numberRows_,&i1,
     498        NULL,permute_,permuteInverse_,array,&numberRows_,&i1,
    516499        NULL,&i0,NULL,integerParameters_,doubleParameters_);
    517500#if 1
  • trunk/ClpFactorization.cpp

    r399 r414  
    8080  int numberRows = model->numberRows();
    8181  int numberColumns = model->numberColumns();
     82  if (!numberRows)
     83    return 0;
    8284  // If too many compressions increase area
    8385  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
     
    287289          // See if worth going sparse and when
    288290          if (numberFtranCounts_>100) {
     291            ftranCountInput_= CoinMax(ftranCountInput_,1.0);
    289292            ftranAverageAfterL_ = CoinMax(ftranCountAfterL_/ftranCountInput_,1.0);
    290293            ftranAverageAfterR_ = CoinMax(ftranCountAfterR_/ftranCountAfterL_,1.0);
    291294            ftranAverageAfterU_ = CoinMax(ftranCountAfterU_/ftranCountAfterR_,1.0);
    292             assert (ftranCountInput_&&ftranCountAfterL_&&ftranCountAfterR_);
    293295            if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
    294296              btranAverageAfterU_ = CoinMax(btranCountAfterU_/btranCountInput_,1.0);
     
    529531  regionSparse->checkClear();
    530532#endif
     533  if (!numberRows_)
     534    return 0;
    531535  if (!networkBasis_) {
    532536    collectStatistics_ = true;
     
    569573    regionSparse->checkClear();
    570574#endif
     575  if (!numberRows_)
     576    return 0;
    571577  if (!networkBasis_) {
    572578    collectStatistics_ = true;
     
    605611                          CoinIndexedVector * regionSparse2) const
    606612{
     613  if (!numberRows_)
     614    return 0;
    607615  if (!networkBasis_) {
    608616    collectStatistics_ = true;
  • trunk/ClpInterior.cpp

    r399 r414  
    644644  assert (!dualR_);
    645645  // create arrays if we are doing KKT
    646   if (cholesky_->type()>20) {
     646  if (cholesky_->type()>=20) {
    647647    primalR_ = new double [nTotal];
    648648    CoinZeroN(primalR_,nTotal);
     
    759759{
    760760  // bad if empty
    761   if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
     761  if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
    762762    problemStatus_=emptyProblem();
    763763    return false;
  • trunk/ClpPredictorCorrector.cpp

    r399 r414  
    8484  //bool firstTime=true;
    8585  //firstFactorization(true);
    86   if (cholesky_->order(this)) {
    87     printf("Not enough memory\n");
     86  int returnCode = cholesky_->order(this);
     87  if (returnCode||cholesky_->symbolic()) {
     88    printf("Error return from symbolic - probably not enough memory\n");
    8889    problemStatus_=4;
    8990    //delete all temporary regions
     
    9394      delete matrix_;
    9495      matrix_ = saveMatrix;
     96    }
     97    // Restore quadratic objective if necessary
     98    if (saveObjective) {
     99      delete objective_;
     100      objective_=saveObjective;
    95101    }
    96102    return -1;
     
    784790  checkSolution();
    785791  CoinMemcpyN(reducedCost_,numberColumns_,dj_);
     792  // If quadratic use last solution
    786793  // Restore quadratic objective if necessary
    787794  if (saveObjective) {
    788795    delete objective_;
    789796    objective_=saveObjective;
     797    objectiveValue_=0.5*(primalObjective_+dualObjective_);
    790798  }
    791799  handler_->message(CLP_BARRIER_END,messages_)
     
    14201428    double maximumRHS;
    14211429    double saveMaximum;
    1422     maximumRHS = maximumAbsElement(deltaY_,numberRows_);
     1430    maximumRHS = CoinMax(maximumAbsElement(deltaY_,numberRows_),1.0e-12);
    14231431    saveMaximum = maximumRHS;
    14241432    if (cholesky_->type()<20) {
  • trunk/ClpSimplex.cpp

    r413 r414  
    54355435        <<CoinMessageEol;
    54365436  } else {
     5437#ifndef NDEBUG   
    54375438    int returnCode=internalFactorize(1);
    54385439    assert (!returnCode);
     5440#else
     5441    internalFactorize(1);
     5442#endif
    54395443  }
    54405444  // put back original costs and then check
  • trunk/ClpSimplexNonlinear.cpp

    r401 r414  
    586586  int nextUnflag=10;
    587587  int nextUnflagIteration=numberIterations_+10;
    588   double nullError=-1.0;
     588  // get two arrays
     589  double * array1 = new double[numberRows_+numberColumns_];
     590  double solutionError=-1.0;
    589591  while (problemStatus_==-1) {
    590592    int result;
     
    624626    pivotRow_=-1;
    625627    result = pivotColumn(rowArray_[3],rowArray_[0],
    626                          columnArray_[0],rowArray_[1],pivotMode,nullError);
     628                         columnArray_[0],rowArray_[1],pivotMode,solutionError,
     629                         array1);
    627630    if (result) {
    628631      if (result==3)
     
    722725    }
    723726  }
     727  delete [] array1;
    724728  return returnCode;
    725729}
     
    729733                                      CoinIndexedVector * spare1, CoinIndexedVector * spare2,
    730734                                      int pivotMode2,
    731                                       double & normFlagged,
     735                                      double & normFlagged,double & normUnflagged,
    732736                                      int & numberNonBasic)
    733737{
     
    742746  sequenceIn_=-1;
    743747  normFlagged=0.0;
     748  normUnflagged=1.0;
    744749  double dualTolerance2 = CoinMin(1.0e-8,1.0e-2*dualTolerance_);
     750  double dualTolerance3 = CoinMin(1.0e-3,1.0e3*dualTolerance_);
    745751  if (!numberNonBasic) {
    746752    //if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
    747753    //printf("infeasible\n");
    748754    if (!pivotMode2||pivotMode2>=10) {
     755      normUnflagged=0.0;
    749756      double bestDj=0.0;
    750757      double bestSuper=0.0;
     
    762769            break;
    763770          case atUpperBound:
    764             if (dj_[iSequence]>dualTolerance_)
     771            if (dj_[iSequence]>dualTolerance3)
    765772              normFlagged += dj_[iSequence]*dj_[iSequence];
    766773            break;
    767774          case atLowerBound:
    768             if (dj_[iSequence]<-dualTolerance_)
     775            if (dj_[iSequence]<-dualTolerance3)
    769776              normFlagged += dj_[iSequence]*dj_[iSequence];
    770777            break;
    771778          case isFree:
    772779          case superBasic:
    773             if (fabs(dj_[iSequence])>dualTolerance_)
     780            if (fabs(dj_[iSequence])>dualTolerance3)
    774781              normFlagged += dj_[iSequence]*dj_[iSequence];
    775782            break;
     
    784791        case atUpperBound:
    785792          if (dj_[iSequence]>dualTolerance_) {
     793            if (dj_[iSequence]>dualTolerance3)
     794              normUnflagged += dj_[iSequence]*dj_[iSequence];
    786795            if (pivotMode2<10) {
    787796              array[iSequence]=-dj_[iSequence];
     
    797806        case atLowerBound:
    798807          if (dj_[iSequence]<-dualTolerance_) {
     808            if (dj_[iSequence]<-dualTolerance3)
     809              normUnflagged += dj_[iSequence]*dj_[iSequence];
    799810            if (pivotMode2<10) {
    800811              array[iSequence]=-dj_[iSequence];
     
    814825#ifndef ALLSUPER
    815826          if (fabs(dj_[iSequence])>dualTolerance_) {
     827            if (fabs(dj_[iSequence])>dualTolerance3)
     828              normUnflagged += dj_[iSequence]*dj_[iSequence];
    816829            nSuper++;
    817830            bestSuper = CoinMax(fabs(dj_[iSequence]),bestSuper);
     
    11421155                                 CoinIndexedVector * spare,
    11431156                                 int & pivotMode,
    1144                                  double & nullError)
     1157                                 double & solutionError,
     1158                                 double * dArray)
    11451159{
    11461160  // say not optimal
     
    11571171  int number=longArray->getNumElements();
    11581172  int numberTotal = numberRows_+numberColumns_;
    1159   double * d = new double [numberTotal];
    11601173  int bestSequence=-1;
    11611174  int bestBasicSequence=-1;
     
    11721185  double djNorm=0.0;
    11731186  double normFlagged=0.0;
     1187  double normUnflagged=0.0;
    11741188  int localPivotMode=pivotMode;
    11751189  bool allFinished=false;
    11761190  bool justOne=false;
    11771191  int returnCode=1;
     1192  double currentObj;
     1193  double predictedObj;
     1194  double thetaObj;
     1195  objective_->stepLength(this,solution_,solution_,0.0,
     1196                                               currentObj,predictedObj,thetaObj);
     1197  double saveObj=currentObj;
     1198#define MINTYPE 1
     1199#if MINTYPE ==2
     1200  // try Shanno's method
     1201  // Copy code from ....cpp.q
     1202#endif
    11781203  // big big loop
    11791204  while (!allFinished) {
     
    11881213#define CONJUGATE 2
    11891214#if CONJUGATE == 0
    1190     bool conjugate = false;
     1215    int conjugate = 0;
    11911216#elif CONJUGATE == 1
    1192     bool conjugate = (pivotMode<10) ? true : false;
     1217    int conjugate = (pivotMode<10) ? MINTYPE : 0;
    11931218#elif CONJUGATE == 2
    1194     bool conjugate = true;
     1219    int conjugate = MINTYPE;
    11951220#else
    1196     bool conjugate = true;
     1221    int conjugate = MINTYPE;
    11971222#endif
    11981223   
     
    12341259        startLocalMode=local;
    12351260      directionVector(longArray,spare,rowArray,local,
    1236                       normFlagged,numberNonBasic);
     1261                      normFlagged,normUnflagged,numberNonBasic);
    12371262      {
    12381263        // check null vector
     
    12521277        }
    12531278#if CLP_DEBUG > 0
    1254         if (handler_->logLevel()>3&&largest>1.0e-8)
    1255           printf("largest non null %g on row %d\n",largest,iLargest);
    1256 #endif
    1257         if (nullError<0.0) {
    1258           nullError=largest;
    1259         } else if (largest>max(1.0e-8,1.0e2*nullError)&&
     1279        if ((handler_->logLevel()&32)&&largest>1.0e-8)
     1280          printf("largest rhs error %g on row %d\n",largest,iLargest);
     1281#endif
     1282        if (solutionError<0.0) {
     1283          solutionError=largest;
     1284        } else if (largest>max(1.0e-8,1.0e2*solutionError)&&
    12601285                   factorization_->pivots()) {
    12611286          longArray->clear();
     
    12821307      }
    12831308#if CONJUGATE == 2
    1284       conjugate = (localPivotMode<10) ? true : false;
     1309      conjugate = (localPivotMode<10) ? MINTYPE : 0;
    12851310#endif
    12861311      if (!nPasses) {
    12871312        djNorm0 = CoinMax(djNorm,1.0e-20);
    1288         memcpy(d,work,numberTotal*sizeof(double));
     1313        memcpy(dArray,work,numberTotal*sizeof(double));
    12891314        if (sequenceIn_>=0&&numberNonBasic==1) {
    12901315          // see if simple move
    1291           double currentObj;
    1292           double predictedObj;
    1293           double thetaObj;
    12941316          double objTheta2 =objective_->stepLength(this,solution_,work,1.0e30,
    12951317                                               currentObj,predictedObj,thetaObj);
     
    13791401            pivotRow_ = -1;
    13801402            theta_=0.0;
    1381             delete [] d;
    13821403            return 0;
    13831404          }
     
    13891410          //     djNorm,djNormSave,beta);
    13901411          for (iSequence=0;iSequence<numberTotal;iSequence++)
    1391             d[iSequence] = work[iSequence] + beta *d[iSequence];
     1412            dArray[iSequence] = work[iSequence] + beta *dArray[iSequence];
    13921413        } else {
    13931414          for (iSequence=0;iSequence<numberTotal;iSequence++)
    1394             d[iSequence] = work[iSequence];
     1415            dArray[iSequence] = work[iSequence];
    13951416        }
    13961417      }
     
    15031524          lastSequenceIn=-1;
    15041525        }
    1505         if (!nFlagged2) {
     1526        if (!nFlagged2||(normFlagged+normUnflagged<1.0e-8)) {
    15061527          primalColumnPivot_->setLooksOptimal(true);
    1507           delete [] d;
    15081528          return 0;
    15091529        } else {
     
    15291549        //if (flagged(iSequence)
    15301550        //  continue;
    1531         double alpha = d[iSequence];
     1551        double alpha = dArray[iSequence];
    15321552        Status thisStatus = getStatus(iSequence);
    15331553        double oldValue = solution_[iSequence];
     
    15761596      theta_ = CoinMin(theta,basicTheta);
    15771597      // Now find minimum of function
    1578       double currentObj;
    1579       double predictedObj;
    1580       double thetaObj;
    1581       double objTheta2 =objective_->stepLength(this,solution_,d,CoinMin(theta,basicTheta),
     1598      double objTheta2 =objective_->stepLength(this,solution_,dArray,CoinMin(theta,basicTheta),
    15821599                                               currentObj,predictedObj,thetaObj);
    15831600#ifdef CLP_DEBUG
     
    15881605        double offset;
    15891606        const double * gradient = objective_->gradient(this,
    1590                                                        d, offset,
     1607                                                       dArray, offset,
    15911608                                                       true,0);
    15921609        double product=0.0;
    15931610        for (iSequence = 0;iSequence<numberColumns_;iSequence++) {
    1594           double alpha = d[iSequence];
     1611          double alpha = dArray[iSequence];
    15951612          double value = alpha*gradient[iSequence];
    15961613          product += value;
     
    15991616#ifdef INCLUDESLACK
    16001617        for (;iSequence<numberColumns_+numberRows_;iSequence++) {
    1601           double alpha = d[iSequence];
     1618          double alpha = dArray[iSequence];
    16021619          double value = alpha*cost_[iSequence];
    16031620          product += value;
     
    16511668      }
    16521669#endif
     1670      if (handler_->logLevel()>3&&objTheta!=theta_) {
     1671        printf("%d pass obj %g,",nPasses,currentObj);
     1672        if (sequenceIn_>=0)
     1673          printf (" Sin %d,",sequenceIn_);
     1674        if (basicTheta==theta_)
     1675          printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
     1676        else
     1677          printf(" basicTheta %g",basicTheta);
     1678        if (theta==theta_)
     1679          printf(" X(%d) non-basicTheta %g",bestSequence,theta);
     1680        else
     1681          printf(" non-basicTheta %g",theta);
     1682        printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
     1683        printf(" djNorm %g\n",djNorm);
     1684      }
    16531685      if (objTheta!=theta_) {
    16541686        //printf("hit boundary after %d passes\n",nPasses);
     
    16651697       
    16661698        //int iSequence = which[iIndex];
    1667         double alpha = d[iSequence];
     1699        double alpha = dArray[iSequence];
    16681700        if (alpha) {
    16691701          double value = solution_[iSequence]+theta_*alpha;
     
    17041736#ifdef CLP_DEBUG
    17051737      if (handler_->logLevel()&32) {
    1706         double currentObj;
    1707         double predictedObj;
    1708         double thetaObj;
    17091738        objective_->stepLength(this,solution_,work,0.0,
    17101739                               currentObj,predictedObj,thetaObj);
     
    17641793        else if (factorization_->pivots()>5)
    17651794          acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
     1795        // should be dArray but seems better this way!
    17661796        double direction = work[bestBasicSequence]>0.0 ? -1.0 : 1.0;
    17671797        // create as packed
     
    18141844              //                  upper_[iSequence]-solution_[iSequence]);
    18151845              double alpha=work2[i];
     1846              // should be dArray but seems better this way!
    18161847              double change=work[iSequence];
    18171848              Status thisStatus = getStatus(iSequence);
     
    19952026    }
    19962027  }
    1997   delete [] d;
    19982028  if (djNorm0<1.0e-12*normFlagged) {
    19992029#ifdef CLP_DEBUG
     
    20022032#endif
    20032033    unflag();
     2034  }
     2035  if (saveObj-currentObj<1.0e-5&&nTotalPasses>2000) {
     2036    normUnflagged=0.0;
     2037    double dualTolerance3 = CoinMin(1.0e-3,1.0e3*dualTolerance_);
     2038    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
     2039      switch(getStatus(iSequence)) {
     2040       
     2041      case basic:
     2042      case ClpSimplex::isFixed:
     2043        break;
     2044      case atUpperBound:
     2045        if (dj_[iSequence]>dualTolerance3)
     2046          normFlagged += dj_[iSequence]*dj_[iSequence];
     2047        break;
     2048      case atLowerBound:
     2049        if (dj_[iSequence]<-dualTolerance3)
     2050          normFlagged += dj_[iSequence]*dj_[iSequence];
     2051        break;
     2052      case isFree:
     2053      case superBasic:
     2054        if (fabs(dj_[iSequence])>dualTolerance3)
     2055          normFlagged += dj_[iSequence]*dj_[iSequence];
     2056        break;
     2057      }
     2058    }
     2059    if (handler_->logLevel()>2)
     2060      printf("possible optimal  %d %d %g %g\n",
     2061             nBigPasses,nTotalPasses,saveObj-currentObj,normFlagged);
     2062    if (normFlagged<1.0e-5) {
     2063      unflag();
     2064      primalColumnPivot_->setLooksOptimal(true);
     2065      returnCode=0;
     2066    }
    20042067  }
    20052068  return returnCode;
     
    23242387  memset(saveRowSolution,0,numberRows*sizeof(double));
    23252388  double * savePi = new double [numberRows];
     2389  double * safeSolution = new double [numberColumns];
    23262390  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
    23272391#define MULTIPLE 0
     
    23992463        }
    24002464      }
     2465      // make sure both accurate
     2466      memset(rowActivity_,0,numberRows_*sizeof(double));
     2467      times(1.0,solution,rowActivity_);
     2468      memset(saveRowSolution,0,numberRows_*sizeof(double));
     2469      times(1.0,saveSolution,saveRowSolution);
    24012470      for (int iRow=0;iRow<numberRows;iRow++) {
    24022471        double alpha =rowActivity_[iRow]-saveRowSolution[iRow];
     
    24442513          solution[iColumn] = lambda * saveSolution[iColumn]
    24452514            + theta * solution[iColumn];
     2515        memset(rowActivity_,0,numberRows_*sizeof(double));
     2516        times(1.0,solution,rowActivity_);
    24462517        if (lambda>0.999) {
    24472518          memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
     
    24532524        int saveScaling = scalingFlag();
    24542525        scaling(0);
     2526        int savePerturbation=perturbation_;
     2527        perturbation_=100;
    24552528        if (saveLogLevel==1)
    24562529          setLogLevel(0);
    2457         int status=startup(0);
     2530        int status=startup(1);
    24582531        if (!status) {
    24592532          int numberTotal = numberRows_+numberColumns_;
     
    25102583            longArray->clear();
    25112584            double normFlagged=0.0;
     2585            double normUnflagged=0.0;
    25122586            int numberNonBasic=0;
    25132587            directionVector(longArray,spare,rowArray,0,
    2514                             normFlagged,numberNonBasic);
     2588                            normFlagged,normUnflagged,numberNonBasic);
     2589            if (normFlagged+normUnflagged<1.0e-8)
     2590              break;  //looks optimal
    25152591            double djNorm=0.0;
    25162592            int iIndex;
     
    26262702        finish();
    26272703        scaling(saveScaling);
     2704        perturbation_=savePerturbation;
    26282705        setLogLevel(saveLogLevel);
    26292706#endif
     2707        // redo rowActivity
     2708        memset(rowActivity_,0,numberRows_*sizeof(double));
     2709        times(1.0,solution,rowActivity_);
    26302710        if (theta>0.99999&&theta2<1.9) {
    26312711          // If big changes then tighten
     
    28362916      ClpSimplex::primal(1);
    28372917      setLogLevel(saveLogLevel);
     2918#ifdef CLP_DEBUG
    28382919      if (this->status()) {
    28392920        writeMps("xx.mps");
    28402921      }
    2841       assert (this->status()!=1); // must be feasible
     2922#endif
     2923      if (this->status()==1) {
     2924        // not feasible ! - backtrack and exit
     2925        // use safe solution
     2926        memcpy(solution,safeSolution,numberColumns*sizeof(double));
     2927        memcpy(saveSolution,solution,numberColumns*sizeof(double));
     2928        memset(rowActivity_,0,numberRows_*sizeof(double));
     2929        times(1.0,solution,rowActivity_);
     2930        memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
     2931        memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
     2932        memcpy(status_,saveStatus,numberRows+numberColumns);
     2933        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     2934          iColumn=listNonLinearColumn[jNon];
     2935          columnLower[iColumn]=CoinMax(solution[iColumn]
     2936                                       -trust[jNon],
     2937                                       trueLower[jNon]);
     2938          columnUpper[iColumn]=CoinMin(solution[iColumn]
     2939                                       +trust[jNon],
     2940                                       trueUpper[jNon]);
     2941        }
     2942        break;
     2943      } else {
     2944        // save in case problems
     2945        memcpy(safeSolution,solution,numberColumns*sizeof(double));
     2946      }
    28422947      if (problemStatus_==3)
    28432948        break; // probably user interrupt
     
    28912996  }
    28922997  delete [] saveSolution;
     2998  delete [] safeSolution;
    28932999  delete [] saveRowSolution;
    28943000  for (iPass=0;iPass<3;iPass++)
  • trunk/ClpSolve.cpp

    r399 r414  
    1616#include "ClpSimplex.hpp"
    1717#include "ClpInterior.hpp"
     18#include "ClpCholeskyDense.hpp"
     19#include "ClpCholeskyBase.hpp"
    1820#ifdef REAL_BARRIER
    1921#include "ClpCholeskyWssmp.hpp"
    2022#include "ClpCholeskyWssmpKKT.hpp"
    2123//#include "ClpCholeskyTaucs.hpp"
     24#include "ClpCholeskyUfl.hpp"
     25#else
     26static int numberBarrier=0;
    2227#endif
    2328#include "ClpSolve.hpp"
     
    102107  int doSprint=0;
    103108  int doSlp=0;
    104   if (method!=ClpSolve::useDual&&method!=ClpSolve::useBarrier) {
     109  if (method!=ClpSolve::useDual&&method!=ClpSolve::useBarrier
     110      &&method!=ClpSolve::useBarrierNoCross) {
    105111    switch (options.getSpecialOption(1)) {
    106112    case 0:
     
    905911    timeX=time2;
    906912    model2->setNumberIterations(model2->numberIterations()+totalIterations);
    907   } else if (method==ClpSolve::useBarrier) {
     913  } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
    908914    //printf("***** experimental pretty crude barrier\n");
    909915    //#define SAVEIT 1
     
    922928    bool aggressiveGamma=false;
    923929    bool scale=false;
     930    bool doKKT=false;
     931    if (barrierOptions&32) {
     932      barrierOptions &= ~32;
     933      doKKT=true;
     934    }
     935    if (barrierOptions&16) {
     936      barrierOptions &= ~16;
     937      aggressiveGamma=true;
     938    }
    924939    if (barrierOptions&8) {
    925940      barrierOptions &= ~8;
    926       aggressiveGamma=true;
    927     }
    928     if (barrierOptions&4) {
    929       barrierOptions &= ~4;
    930941      scale=true;
    931942    }
     
    940951#ifdef REAL_BARRIER
    941952    if (quadraticObj) {
    942       barrierOptions = 2;
     953      doKKT=true;
    943954    }
    944955    switch (barrierOptions) {
    945956    case 0:
    946       {
    947         ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
     957      if (!doKKT) {
     958        ClpCholeskyBase * cholesky = new ClpCholeskyBase(100);
    948959        barrier.setCholesky(cholesky);
     960      } else {
     961        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
     962        cholesky->setKKT(true);
     963        barrier.setCholesky(cholesky);
    949964      }
    950965      break;
    951966    case 1:
     967      if (!doKKT) {
     968        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
     969        barrier.setCholesky(cholesky);
     970      } else {
     971        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
     972        cholesky->setKKT(true);
     973        barrier.setCholesky(cholesky);
     974      }
     975      break;
     976    case 2:
    952977      {
    953978        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
    954979        barrier.setCholesky(cholesky);
    955       }
    956       break;
    957     case 2:
    958       {
     980        assert (!doKKT);
     981      }
     982      break;
     983    case 3:
     984      if (!doKKT) {
     985        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
     986        barrier.setCholesky(cholesky);
     987      } else {
    959988        ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
    960989        barrier.setCholesky(cholesky);
    961990      }
    962991      break;
    963     case 3:
     992    case 4:
     993      if (!doKKT) {
     994        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl(100);
     995        barrier.setCholesky(cholesky);
     996      } else {
     997        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
     998        cholesky->setKKT(true);
     999        barrier.setCholesky(cholesky);
     1000      }
    9641001      break;
    9651002    default:
     
    9691006    //ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
    9701007    //barrier.setCholesky(cholesky);
     1008#else
     1009    if (quadraticObj) {
     1010      doKKT=true;
     1011    }
     1012    switch (barrierOptions) {
     1013    case 1:
     1014      if (!doKKT) {
     1015        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
     1016        barrier.setCholesky(cholesky);
     1017      } else {
     1018        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
     1019        cholesky->setKKT(true);
     1020        barrier.setCholesky(cholesky);
     1021      }
     1022      break;
     1023    default:
     1024      if (!numberBarrier)
     1025        std::cout<<"Warning - the default ordering is just on row counts! "
     1026                 <<"The factorization is being improved"<<std::endl;
     1027      numberBarrier++;
     1028      if (!doKKT) {
     1029        ClpCholeskyBase * cholesky = new ClpCholeskyBase(100);
     1030        barrier.setCholesky(cholesky);
     1031      } else {
     1032        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
     1033        cholesky->setKKT(true);
     1034        barrier.setCholesky(cholesky);
     1035      }
     1036      break;
     1037    }
    9711038#endif
    9721039    int numberRows = model2->numberRows();
     
    10781145      }
    10791146    }
    1080     if (maxIts&&barrierStatus<4&&!quadraticObj) {
    1081       //printf("***** crossover - needs more thought on difficult models\n");
     1147    if (method==ClpSolve::useBarrier) {
     1148      if (maxIts&&barrierStatus<4&&!quadraticObj) {
     1149        //printf("***** crossover - needs more thought on difficult models\n");
    10821150#if SAVEIT==1
    1083       model2->ClpSimplex::saveModel("xx.save");
     1151        model2->ClpSimplex::saveModel("xx.save");
    10841152#endif
    1085       // make sure no status left
    1086       model2->createStatus();
    1087       // solve
    1088       model2->setPerturbation(100);
     1153        // make sure no status left
     1154        model2->createStatus();
     1155        // solve
     1156        model2->setPerturbation(100);
    10891157#if 1
    1090       // throw some into basis
    1091       {
    1092         int numberRows = model2->numberRows();
    1093         int numberColumns = model2->numberColumns();
    1094         double * dsort = new double[numberColumns];
    1095         int * sort = new int[numberColumns];
    1096         int n=0;
    1097         const double * columnLower = model2->columnLower();
    1098         const double * columnUpper = model2->columnUpper();
    1099         double * primalSolution = model2->primalColumnSolution();
    1100         const double * dualSolution = model2->dualColumnSolution();
    1101         double tolerance = 10.0*primalTolerance_;
    1102         int i;
    1103         for ( i=0;i<numberRows;i++)
    1104           model2->setRowStatus(i,superBasic);
    1105         for ( i=0;i<numberColumns;i++) {
    1106           double distance = CoinMin(columnUpper[i]-primalSolution[i],
    1107                                 primalSolution[i]-columnLower[i]);
    1108           if (distance>tolerance) {
    1109             if (fabs(dualSolution[i])<1.0e-5)
    1110               distance *= 100.0;
    1111             dsort[n]=-distance;
    1112             sort[n++]=i;
    1113             model2->setStatus(i,superBasic);
    1114           } else if (distance>primalTolerance_) {
    1115             model2->setStatus(i,superBasic);
    1116           } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
    1117             model2->setStatus(i,atLowerBound);
    1118             primalSolution[i]=columnLower[i];
    1119           } else {
    1120             model2->setStatus(i,atUpperBound);
    1121             primalSolution[i]=columnUpper[i];
     1158        // throw some into basis
     1159        {
     1160          int numberRows = model2->numberRows();
     1161          int numberColumns = model2->numberColumns();
     1162          double * dsort = new double[numberColumns];
     1163          int * sort = new int[numberColumns];
     1164          int n=0;
     1165          const double * columnLower = model2->columnLower();
     1166          const double * columnUpper = model2->columnUpper();
     1167          double * primalSolution = model2->primalColumnSolution();
     1168          const double * dualSolution = model2->dualColumnSolution();
     1169          double tolerance = 10.0*primalTolerance_;
     1170          int i;
     1171          for ( i=0;i<numberRows;i++)
     1172            model2->setRowStatus(i,superBasic);
     1173          for ( i=0;i<numberColumns;i++) {
     1174            double distance = CoinMin(columnUpper[i]-primalSolution[i],
     1175                                      primalSolution[i]-columnLower[i]);
     1176            if (distance>tolerance) {
     1177              if (fabs(dualSolution[i])<1.0e-5)
     1178                distance *= 100.0;
     1179              dsort[n]=-distance;
     1180              sort[n++]=i;
     1181              model2->setStatus(i,superBasic);
     1182            } else if (distance>primalTolerance_) {
     1183              model2->setStatus(i,superBasic);
     1184            } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
     1185              model2->setStatus(i,atLowerBound);
     1186              primalSolution[i]=columnLower[i];
     1187            } else {
     1188              model2->setStatus(i,atUpperBound);
     1189              primalSolution[i]=columnUpper[i];
     1190            }
    11221191          }
     1192          CoinSort_2(dsort,dsort+n,sort);
     1193          n = CoinMin(numberRows,n);
     1194          for ( i=0;i<n;i++) {
     1195            int iColumn = sort[i];
     1196            model2->setStatus(iColumn,basic);
     1197          }
     1198          delete [] sort;
     1199          delete [] dsort;
    11231200        }
    1124         CoinSort_2(dsort,dsort+n,sort);
    1125         n = CoinMin(numberRows,n);
    1126         for ( i=0;i<n;i++) {
    1127           int iColumn = sort[i];
    1128           model2->setStatus(iColumn,basic);
     1201        if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
     1202          int numberRows = model2->numberRows();
     1203          int numberColumns = model2->numberColumns();
     1204          // just primal values pass
     1205          double saveScale = model2->objectiveScale();
     1206          model2->setObjectiveScale(1.0e-3);
     1207          model2->primal(2);
     1208          model2->setObjectiveScale(saveScale);
     1209          // save primal solution and copy back dual
     1210          CoinMemcpyN(model2->primalRowSolution(),
     1211                      numberRows,rowPrimal);
     1212          CoinMemcpyN(rowDual,
     1213                      numberRows,model2->dualRowSolution());
     1214          CoinMemcpyN(model2->primalColumnSolution(),
     1215                      numberColumns,columnPrimal);
     1216          CoinMemcpyN(columnDual,
     1217                      numberColumns,model2->dualColumnSolution());
     1218          //model2->primal(1);
     1219          // clean up reduced costs and flag variables
     1220          {
     1221            double * dj = model2->dualColumnSolution();
     1222            double * cost = model2->objective();
     1223            double * saveCost = new double[numberColumns];
     1224            memcpy(saveCost,cost,numberColumns*sizeof(double));
     1225            double * saveLower = new double[numberColumns];
     1226            double * lower = model2->columnLower();
     1227            memcpy(saveLower,lower,numberColumns*sizeof(double));
     1228            double * saveUpper = new double[numberColumns];
     1229            double * upper = model2->columnUpper();
     1230            memcpy(saveUpper,upper,numberColumns*sizeof(double));
     1231            int i;
     1232            double tolerance = 10.0*dualTolerance_;
     1233            for ( i=0;i<numberColumns;i++) {
     1234              if (model2->getStatus(i)==basic) {
     1235                dj[i]=0.0;
     1236              } else if (model2->getStatus(i)==atLowerBound) {
     1237                if (optimizationDirection_*dj[i]<tolerance) {
     1238                  if (optimizationDirection_*dj[i]<0.0) {
     1239                    //if (dj[i]<-1.0e-3)
     1240                    //printf("bad dj at lb %d %g\n",i,dj[i]);
     1241                    cost[i] -= dj[i];
     1242                    dj[i]=0.0;
     1243                  }
     1244                } else {
     1245                  upper[i]=lower[i];
     1246                }
     1247              } else if (model2->getStatus(i)==atUpperBound) {
     1248                if (optimizationDirection_*dj[i]>tolerance) {
     1249                  if (optimizationDirection_*dj[i]>0.0) {
     1250                    //if (dj[i]>1.0e-3)
     1251                    //printf("bad dj at ub %d %g\n",i,dj[i]);
     1252                    cost[i] -= dj[i];
     1253                    dj[i]=0.0;
     1254                  }
     1255                } else {
     1256                  lower[i]=upper[i];
     1257                }
     1258              }
     1259            }
     1260            // just dual values pass
     1261            //model2->setLogLevel(63);
     1262            //model2->setFactorizationFrequency(1);
     1263            model2->dual(2);
     1264            memcpy(cost,saveCost,numberColumns*sizeof(double));
     1265            delete [] saveCost;
     1266            memcpy(lower,saveLower,numberColumns*sizeof(double));
     1267            delete [] saveLower;
     1268            memcpy(upper,saveUpper,numberColumns*sizeof(double));
     1269            delete [] saveUpper;
     1270          }
     1271          // and finish
     1272          // move solutions
     1273          CoinMemcpyN(rowPrimal,
     1274                      numberRows,model2->primalRowSolution());
     1275          CoinMemcpyN(columnPrimal,
     1276                      numberColumns,model2->primalColumnSolution());
    11291277        }
    1130         delete [] sort;
    1131         delete [] dsort;
    1132       }
    1133       if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
    1134         int numberRows = model2->numberRows();
    1135         int numberColumns = model2->numberColumns();
    1136         // just primal values pass
    11371278        double saveScale = model2->objectiveScale();
    11381279        model2->setObjectiveScale(1.0e-3);
    11391280        model2->primal(2);
    11401281        model2->setObjectiveScale(saveScale);
    1141         // save primal solution and copy back dual
    1142         CoinMemcpyN(model2->primalRowSolution(),
    1143                     numberRows,rowPrimal);
    1144         CoinMemcpyN(rowDual,
    1145                     numberRows,model2->dualRowSolution());
    1146         CoinMemcpyN(model2->primalColumnSolution(),
    1147                     numberColumns,columnPrimal);
    1148         CoinMemcpyN(columnDual,
    1149                     numberColumns,model2->dualColumnSolution());
    1150         //model2->primal(1);
    1151         // clean up reduced costs and flag variables
    1152         {
    1153           double * dj = model2->dualColumnSolution();
    1154           double * cost = model2->objective();
    1155           double * saveCost = new double[numberColumns];
    1156           memcpy(saveCost,cost,numberColumns*sizeof(double));
    1157           double * saveLower = new double[numberColumns];
    1158           double * lower = model2->columnLower();
    1159           memcpy(saveLower,lower,numberColumns*sizeof(double));
    1160           double * saveUpper = new double[numberColumns];
    1161           double * upper = model2->columnUpper();
    1162           memcpy(saveUpper,upper,numberColumns*sizeof(double));
    1163           int i;
    1164           double tolerance = 10.0*dualTolerance_;
    1165           for ( i=0;i<numberColumns;i++) {
    1166             if (model2->getStatus(i)==basic) {
    1167               dj[i]=0.0;
    1168             } else if (model2->getStatus(i)==atLowerBound) {
    1169               if (optimizationDirection_*dj[i]<tolerance) {
    1170                 if (optimizationDirection_*dj[i]<0.0) {
    1171                   //if (dj[i]<-1.0e-3)
    1172                   //printf("bad dj at lb %d %g\n",i,dj[i]);
    1173                   cost[i] -= dj[i];
    1174                   dj[i]=0.0;
    1175                 }
    1176               } else {
    1177                 upper[i]=lower[i];
    1178               }
    1179             } else if (model2->getStatus(i)==atUpperBound) {
    1180               if (optimizationDirection_*dj[i]>tolerance) {
    1181                 if (optimizationDirection_*dj[i]>0.0) {
    1182                   //if (dj[i]>1.0e-3)
    1183                   //printf("bad dj at ub %d %g\n",i,dj[i]);
    1184                   cost[i] -= dj[i];
    1185                   dj[i]=0.0;
    1186                 }
    1187               } else {
    1188                 lower[i]=upper[i];
    1189               }
    1190             }
    1191           }
    1192           // just dual values pass
    1193           //model2->setLogLevel(63);
    1194           //model2->setFactorizationFrequency(1);
    1195           model2->dual(2);
    1196           memcpy(cost,saveCost,numberColumns*sizeof(double));
    1197           delete [] saveCost;
    1198           memcpy(lower,saveLower,numberColumns*sizeof(double));
    1199           delete [] saveLower;
    1200           memcpy(upper,saveUpper,numberColumns*sizeof(double));
    1201           delete [] saveUpper;
    1202         }
    1203         // and finish
    1204         // move solutions
    1205         CoinMemcpyN(rowPrimal,
    1206                     numberRows,model2->primalRowSolution());
    1207         CoinMemcpyN(columnPrimal,
    1208                     numberColumns,model2->primalColumnSolution());
    1209       }
    1210       double saveScale = model2->objectiveScale();
    1211       model2->setObjectiveScale(1.0e-3);
    1212       model2->primal(2);
    1213       model2->setObjectiveScale(saveScale);
    1214       model2->primal(1);
     1282        model2->primal(1);
    12151283#else
    1216       // just primal
    1217       model2->primal(1);
     1284        // just primal
     1285        model2->primal(1);
    12181286#endif
    1219     } else if (barrierStatus==4) {
    1220       // memory problems
    1221       model2->setPerturbation(savePerturbation);
    1222       model2->createStatus();
    1223       model2->dual();
    1224     } else if (maxIts&&quadraticObj) {
    1225       // make sure no status left
    1226       model2->createStatus();
    1227       // solve
    1228       model2->setPerturbation(100);
    1229       model2->reducedGradient(1);
     1287      } else if (barrierStatus==4) {
     1288        // memory problems
     1289        model2->setPerturbation(savePerturbation);
     1290        model2->createStatus();
     1291        model2->dual();
     1292      } else if (maxIts&&quadraticObj) {
     1293        // make sure no status left
     1294        model2->createStatus();
     1295        // solve
     1296        model2->setPerturbation(100);
     1297        model2->reducedGradient(1);
     1298      }
    12301299    }
    12311300    model2->setMaximumIterations(saveMaxIts);
  • trunk/Test/ClpMain.cpp

    r398 r414  
    1414#include "CoinPragma.hpp"
    1515#include "CoinHelperFunctions.hpp"
    16 #define CLPVERSION "0.99.8"
     16#define CLPVERSION "0.99.9"
    1717
    1818#include "CoinMpsIO.hpp"
     
    6363  DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    6464  PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE,
    65   CHOLESKY,BARRIERSCALE,GAMMA,
     65  CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,
    6666 
    6767  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,
     
    494494        length1 = thisOne.length();
    495495      }
    496       if (check.length()<=length1) {
     496      if (check.length()<=length1&&length2<=check.length()) {
    497497        unsigned int i;
    498498        for (i=0;i<check.length();i++) {
     
    927927      (
    928928       "This command solves the current model using the  primal dual predictor \
    929 corrector algorithm.  This is not a sophisticated version just something JJF \
    930 knocked up"
     929corrector algorithm."
    931930
    932931       );
     
    943942    parameters[numberParameters++]=
    944943      ClpItem("chol!esky","Which cholesky algorithm",
    945               "wssmp",CHOLESKY,false);
     944              "base",CHOLESKY,false);
     945    parameters[numberParameters-1].append("dense");
     946#ifdef REAL_BARRIER
    946947    parameters[numberParameters-1].append("fudge!Long");
    947     parameters[numberParameters-1].append("KKT");
    948     parameters[numberParameters-1].append("dense");
     948    parameters[numberParameters-1].append("wssmp");
     949    parameters[numberParameters-1].append("Uni!versityOfFlorida");
     950#endif
    949951    parameters[numberParameters++]=
    950952      ClpItem("crash","Whether to create basis for problem",
     
    956958 variables into basis with the aim of getting dual feasible.  On the whole dual seems to be\
    957959 better without it and there alernative types of 'crash' for primal e.g. 'idiot' or 'sprint'."
     960       );
     961    parameters[numberParameters++]=
     962      ClpItem("cross!over","Whether to get a basic solution after barrier",
     963              "on",CROSSOVER);
     964    parameters[numberParameters-1].append("off");
     965    parameters[numberParameters-1].append("maybe");
     966    parameters[numberParameters-1].setLonghelp
     967      (
     968       "Interior point algorithms do not obtain a basic solution (and \
     969the feasibility criterion is a bit suspect (JJF)).  This option will crossover \
     970to a basic solution suitable for ranging or branch and cut.  With the current state \
     971of quadratic it may be a good idea to switch off crossover for quadratic (and maybe \
     972presolve as well) - the option maybe does this."
    958973       );
    959974    parameters[numberParameters++]=
     
    11021117       );
    11031118    parameters[numberParameters++]=
     1119      ClpItem("KKT","Whether to use KKT factorization",
     1120              "off",KKT,false);
     1121    parameters[numberParameters-1].append("on");
     1122    parameters[numberParameters++]=
    11041123      ClpItem("log!Level","Level of detail in output",
    11051124              -1,63,LOGLEVEL);
     
    11241143    parameters[numberParameters-1].setLonghelp
    11251144      (
    1126        "If this is at its default value of 200 then in this executable clp will guess at a\
     1145       "If this is at its initial value of 201 then in this executable clp will guess at a\
    11271146 value to use.  Otherwise the user can set a value.  The code may decide to re-factorize\
    11281147 earlier."
     
    14441463    int gamma=0;
    14451464    int scaleBarrier=0;
     1465    int doKKT=0;
     1466    int crossover=2; // do crossover unless quadratic
    14461467   
    14471468    int iModel=0;
     
    17301751              scaleBarrier=action;
    17311752              break;
     1753            case KKT:
     1754              doKKT=action;
     1755              break;
    17321756            default:
    17331757              abort();
     
    17541778                presolveType=ClpSolve::presolveOff;
    17551779              solveOptions.setPresolveType(presolveType,preSolve);
    1756               if (type==DUALSIMPLEX)
     1780              if (type==DUALSIMPLEX) {
    17571781                method=ClpSolve::useDual;
    1758               else if (type==PRIMALSIMPLEX)
     1782              } else if (type==PRIMALSIMPLEX) {
    17591783                method=ClpSolve::usePrimalorSprint;
    1760               else
     1784              } else {
    17611785                method = ClpSolve::useBarrier;
     1786                if (crossover==1) {
     1787                  method=ClpSolve::useBarrierNoCross;
     1788                } else if (crossover==2) {
     1789                  ClpObjective * obj = models[iModel].objectiveAsObject();
     1790                  if (obj->type()>1) {
     1791                    method=ClpSolve::useBarrierNoCross;
     1792                    presolveType=ClpSolve::presolveOff;
     1793                    solveOptions.setPresolveType(presolveType,preSolve);
     1794                  }
     1795                }
     1796              }
    17621797              solveOptions.setSolveType(method);
    17631798              if (method==ClpSolve::useDual) {
     
    17971832                  }
    17981833                }
    1799               } else if (method==ClpSolve::useBarrier) {
     1834              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
    18001835                int barrierOptions = choleskyType;
    18011836                if (scaleBarrier)
    1802                   barrierOptions |= 4;
     1837                  barrierOptions |= 8;
    18031838                if (gamma)
    1804                   barrierOptions |= 8;
     1839                  barrierOptions |= 16;
     1840                if (doKKT)
     1841                  barrierOptions |= 32;
    18051842                solveOptions.setSpecialOption(4,barrierOptions);
    18061843              }
  • trunk/include/ClpCholeskyBase.hpp

    r389 r414  
    55
    66#include "CoinPragma.hpp"
     7#include "CoinFinite.hpp"
    78
    89class ClpInterior;
    910
    10 /** Abstract Base class for Clp Cholesky factorization
     11/** Base class for Clp Cholesky factorization
     12    Will do better factorization.  very crude ordering
    1113
    1214    Derived classes will be using sophisticated methods apart from ClpCholeskyDense
    1315*/
    1416
     17class ClpCholeskyDense;
     18class ClpMatrixBase;
    1519class ClpCholeskyBase  {
    1620 
    1721public:
    18    /**@name Virtual methods that the derived classes must provide  */
     22   /**@name Virtual methods that the derived classes may provide  */
    1923   //@{
    2024  /** Orders rows and saves pointer to matrix.and model.
    21    returns non-zero if not enough memory */
    22   virtual int order(ClpInterior * model) = 0;
     25   returns non-zero if not enough memory.
     26   You can use preOrder to set up ADAT
     27   If using default symbolic etc then must set sizeFactor_ to
     28   size of input matrix to order (and to symbolic).
     29   Also just permute_ and permuteInverse_ should be created */
     30  virtual int order(ClpInterior * model);
     31  /** Does Symbolic factorization given permutation.
     32      This is called immediately after order.  If user provides this then
     33      user must provide factorize and solve.  Otherwise the default factorization is used
     34      returns non-zero if not enough memory */
     35  virtual int symbolic();
    2336  /** Factorize - filling in rowsDropped and returning number dropped.
    2437      If return code negative then out of memory */
    25   virtual int factorize(const double * diagonal, int * rowsDropped) =0;
     38  virtual int factorize(const double * diagonal, int * rowsDropped) ;
    2639  /** Uses factorization to solve. */
    27   virtual void solve (double * region) = 0;
     40  virtual void solve (double * region) ;
    2841  /** Uses factorization to solve. - given as if KKT.
    2942   region1 is rows+columns, region2 is rows */
     
    5467  inline int numberRows() const
    5568  {return numberRows_;};
     69  /// Return size
     70  inline CoinBigIndex size() const
     71  { return sizeFactor_;};
     72  /// Return sparseFactor
     73  inline double * sparseFactor() const
     74  { return sparseFactor_;};
     75  /// Return diagonal
     76  inline double * diagonal() const
     77  { return diagonal_;};
     78  /// Return workDouble
     79  inline double * workDouble() const
     80  { return workDouble_;};
     81  /// If KKT on
     82  inline bool kkt() const
     83  { return doKKT_;};
     84  /// Set KKT
     85  inline void setKKT(bool yesNo)
     86  { doKKT_ = yesNo;};
    5687   //@}
    5788 
    5889 
    59 protected:
    60 
    61    /**@name Constructors, destructor<br>
    62       <strong>NOTE</strong>: All constructors are protected. There's no need
    63       to expose them, after all, this is an abstract class. */
     90public:
     91
     92   /**@name Constructors, destructor
     93    */
    6494   //@{
    65    /** Default constructor. */
    66    ClpCholeskyBase();
     95  /** Constructor which has dense columns activated.
     96      Default is off. */
     97  ClpCholeskyBase(int denseThreshold=-1);
    6798   /** Destructor (has to be public) */
    68 public:
    6999   virtual ~ClpCholeskyBase();
    70 protected:
    71100  // Copy
    72101   ClpCholeskyBase(const ClpCholeskyBase&);
     
    77106  ///@name Other
    78107  /// Clone
    79 public:
    80   virtual ClpCholeskyBase * clone() const = 0;
     108  virtual ClpCholeskyBase * clone() const;
    81109 
    82110  /// Returns type
    83111  inline int type() const
    84   { return type_;};
     112  { if (doKKT_) return 100; else return type_;};
    85113protected:
    86114  /// Sets type
     
    88116   //@}
    89117   
     118  /**@name Symbolic, factor and solve */
     119  //@{
     120  /** Symbolic1  - works out size without clever stuff.
     121      Uses upper triangular as much easier.
     122      Returns size
     123   */
     124  int symbolic1(const CoinBigIndex * Astart, const int * Arow);
     125  /** Symbolic2  - Fills in indices
     126      Uses lower triangular so can do cliques etc
     127   */
     128  void symbolic2(const CoinBigIndex * Astart, const int * Arow);
     129  /** Factorize - filling in rowsDropped and returning number dropped
     130      in integerParam.
     131   */
     132  void factorizePart2(int * rowsDropped) ;
     133  /** solve - 1 just first half, 2 just second half - 3 both.
     134  If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
     135  */
     136  void solve(double * region, int type);
     137  /// Forms ADAT - returns nonzero if not enough memory
     138  int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT);
     139  //@}
    90140   
    91141protected:
    92    /**@name Data members
    93       The data members are protected to allow access for derived classes. */
    94    //@{
    95    /// type (may be useful) if > 20 do KKT
    96    int type_;
    97   /// pivotTolerance.
    98   double pivotTolerance_;
     142  /**@name Data members
     143     The data members are protected to allow access for derived classes. */
     144  //@{
     145  /// type (may be useful) if > 20 do KKT
     146  int type_;
     147  /// Doing full KKT (only used if default symbolic and factorization)
     148  bool doKKT_;
    99149  /// zeroTolerance.
    100150  double zeroTolerance_;
     
    111161  /// rowsDropped
    112162  char * rowsDropped_;
    113   /// permuteIn.
    114   int * permuteIn_;
    115   /// permuteOut.
    116   int * permuteOut_;
     163  /// permute inverse.
     164  int * permuteInverse_;
     165  /// main permute.
     166  int * permute_;
    117167  /// numberRowsDropped.  Number of rows gone
    118168  int numberRowsDropped_;
     169  /// sparseFactor.
     170  double * sparseFactor_;
     171  /// choleskyStart - element starts
     172  CoinBigIndex * choleskyStart_;
     173  /// choleskyRow (can be shorter than sparsefactor)
     174  int * choleskyRow_;
     175  /// Index starts
     176  CoinBigIndex * indexStart_;
     177  /// Diagonal
     178  double * diagonal_;
     179  /// double work array
     180  double * workDouble_;
     181  /// link array
     182  int * link_;
     183  // Integer work array
     184  CoinBigIndex * workInteger_;
     185  // Clique information
     186  int * clique_;
     187  /// sizeFactor.
     188  CoinBigIndex sizeFactor_;
     189  /// Size of index array
     190  CoinBigIndex sizeIndex_;
     191  /// First dense row
     192  int firstDense_;
     193  /// integerParameters
     194  int integerParameters_[64];
     195  /// doubleParameters;
     196  double doubleParameters_[64];
     197  /// Row copy of matrix
     198  ClpMatrixBase * rowCopy_;
     199  /// Dense indicators
     200  char * whichDense_;
     201  /// Dense columns (updated)
     202  double * denseColumn_;
     203  /// Dense cholesky
     204  ClpCholeskyDense * dense_;