source: trunk/Clp/src/ClpCholeskyBase.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 3 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 128.3 KB
Line 
1/* $Id: ClpCholeskyBase.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6/*----------------------------------------------------------------------------*/
7/*      Ordering code - courtesy of Anshul Gupta                              */
8/*      (C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
9/*----------------------------------------------------------------------------*/
10
11/*  A compact no-frills Approximate Minimum Local Fill ordering code.
12
13    References:
14
15[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
16    Edward Rothberg, SGI Manuscript, April 1996.
17[2] An Approximate Minimum Degree Ordering Algorithm.
18    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
19    University of Florida, December 1994.
20*/
21/*----------------------------------------------------------------------------*/
22
23#include "CoinPragma.hpp"
24
25#include <iostream>
26
27#include "ClpCholeskyBase.hpp"
28#include "ClpInterior.hpp"
29#include "ClpHelperFunctions.hpp"
30#include "CoinHelperFunctions.hpp"
31#include "CoinSort.hpp"
32#include "ClpCholeskyDense.hpp"
33#include "ClpMessage.hpp"
34#include "ClpQuadraticObjective.hpp"
35
36//#############################################################################
37// Constructors / Destructor / Assignment
38//#############################################################################
39
40//-------------------------------------------------------------------
41// Default Constructor
42//-------------------------------------------------------------------
43ClpCholeskyBase::ClpCholeskyBase(int denseThreshold)
44  : type_(0)
45  , doKKT_(false)
46  , goDense_(0.7)
47  , choleskyCondition_(0.0)
48  , model_(NULL)
49  , numberTrials_()
50  , numberRows_(0)
51  , status_(0)
52  , rowsDropped_(NULL)
53  , permuteInverse_(NULL)
54  , permute_(NULL)
55  , numberRowsDropped_(0)
56  , sparseFactor_(NULL)
57  , choleskyStart_(NULL)
58  , choleskyRow_(NULL)
59  , indexStart_(NULL)
60  , diagonal_(NULL)
61  , workDouble_(NULL)
62  , link_(NULL)
63  , workInteger_(NULL)
64  , clique_(NULL)
65  , sizeFactor_(0)
66  , sizeIndex_(0)
67  , firstDense_(0)
68  , rowCopy_(NULL)
69  , whichDense_(NULL)
70  , denseColumn_(NULL)
71  , dense_(NULL)
72  , denseThreshold_(denseThreshold)
73{
74  memset(integerParameters_, 0, 64 * sizeof(int));
75  memset(doubleParameters_, 0, 64 * sizeof(double));
76}
77
78//-------------------------------------------------------------------
79// Copy constructor
80//-------------------------------------------------------------------
81ClpCholeskyBase::ClpCholeskyBase(const ClpCholeskyBase &rhs)
82  : type_(rhs.type_)
83  , doKKT_(rhs.doKKT_)
84  , goDense_(rhs.goDense_)
85  , choleskyCondition_(rhs.choleskyCondition_)
86  , model_(rhs.model_)
87  , numberTrials_(rhs.numberTrials_)
88  , numberRows_(rhs.numberRows_)
89  , status_(rhs.status_)
90  , numberRowsDropped_(rhs.numberRowsDropped_)
91{
92  rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
93  permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
94  permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
95  sizeFactor_ = rhs.sizeFactor_;
96  sizeIndex_ = rhs.sizeIndex_;
97  firstDense_ = rhs.firstDense_;
98  sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
99  choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
100  indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
101  choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
102  diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
103#if CLP_LONG_CHOLESKY != 1
104  workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
105#else
106  // actually long double
107  workDouble_ = reinterpret_cast< double * >(ClpCopyOfArray(reinterpret_cast< CoinWorkDouble * >(rhs.workDouble_), numberRows_));
108#endif
109  link_ = ClpCopyOfArray(rhs.link_, numberRows_);
110  workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
111  clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
112  CoinMemcpyN(rhs.integerParameters_, 64, integerParameters_);
113  CoinMemcpyN(rhs.doubleParameters_, 64, doubleParameters_);
114  rowCopy_ = rhs.rowCopy_->clone();
115  whichDense_ = NULL;
116  denseColumn_ = NULL;
117  dense_ = NULL;
118  denseThreshold_ = rhs.denseThreshold_;
119}
120
121//-------------------------------------------------------------------
122// Destructor
123//-------------------------------------------------------------------
124ClpCholeskyBase::~ClpCholeskyBase()
125{
126  delete[] rowsDropped_;
127  delete[] permuteInverse_;
128  delete[] permute_;
129  delete[] sparseFactor_;
130  delete[] choleskyStart_;
131  delete[] choleskyRow_;
132  delete[] indexStart_;
133  delete[] diagonal_;
134  delete[] workDouble_;
135  delete[] link_;
136  delete[] workInteger_;
137  delete[] clique_;
138  delete rowCopy_;
139  delete[] whichDense_;
140  delete[] denseColumn_;
141  delete dense_;
142}
143
144//----------------------------------------------------------------
145// Assignment operator
146//-------------------------------------------------------------------
147ClpCholeskyBase &
148ClpCholeskyBase::operator=(const ClpCholeskyBase &rhs)
149{
150  if (this != &rhs) {
151    type_ = rhs.type_;
152    doKKT_ = rhs.doKKT_;
153    goDense_ = rhs.goDense_;
154    choleskyCondition_ = rhs.choleskyCondition_;
155    model_ = rhs.model_;
156    numberTrials_ = rhs.numberTrials_;
157    numberRows_ = rhs.numberRows_;
158    status_ = rhs.status_;
159    numberRowsDropped_ = rhs.numberRowsDropped_;
160    delete[] rowsDropped_;
161    delete[] permuteInverse_;
162    delete[] permute_;
163    delete[] sparseFactor_;
164    delete[] choleskyStart_;
165    delete[] choleskyRow_;
166    delete[] indexStart_;
167    delete[] diagonal_;
168    delete[] workDouble_;
169    delete[] link_;
170    delete[] workInteger_;
171    delete[] clique_;
172    delete rowCopy_;
173    delete[] whichDense_;
174    delete[] denseColumn_;
175    delete dense_;
176    rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_, numberRows_);
177    permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_, numberRows_);
178    permute_ = ClpCopyOfArray(rhs.permute_, numberRows_);
179    sizeFactor_ = rhs.sizeFactor_;
180    sizeIndex_ = rhs.sizeIndex_;
181    firstDense_ = rhs.firstDense_;
182    sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_, rhs.sizeFactor_);
183    choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_, numberRows_ + 1);
184    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, rhs.sizeFactor_);
185    indexStart_ = ClpCopyOfArray(rhs.indexStart_, numberRows_);
186    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_, sizeIndex_);
187    diagonal_ = ClpCopyOfArray(rhs.diagonal_, numberRows_);
188#if CLP_LONG_CHOLESKY != 1
189    workDouble_ = ClpCopyOfArray(rhs.workDouble_, numberRows_);
190#else
191    // actually long double
192    workDouble_ = reinterpret_cast< double * >(ClpCopyOfArray(reinterpret_cast< CoinWorkDouble * >(rhs.workDouble_), numberRows_));
193#endif
194    link_ = ClpCopyOfArray(rhs.link_, numberRows_);
195    workInteger_ = ClpCopyOfArray(rhs.workInteger_, numberRows_);
196    clique_ = ClpCopyOfArray(rhs.clique_, numberRows_);
197    rowCopy_ = rhs.rowCopy_->clone();
198    whichDense_ = NULL;
199    denseColumn_ = NULL;
200    dense_ = NULL;
201    denseThreshold_ = rhs.denseThreshold_;
202  }
203  return *this;
204}
205// reset numberRowsDropped and rowsDropped.
206void ClpCholeskyBase::resetRowsDropped()
207{
208  numberRowsDropped_ = 0;
209  memset(rowsDropped_, 0, numberRows_);
210}
211/* Uses factorization to solve. - given as if KKT.
212   region1 is rows+columns, region2 is rows */
213void ClpCholeskyBase::solveKKT(CoinWorkDouble *region1, CoinWorkDouble *region2, const CoinWorkDouble *diagonal,
214  CoinWorkDouble diagonalScaleFactor)
215{
216  if (!doKKT_) {
217    int iColumn;
218    int numberColumns = model_->numberColumns();
219    int numberTotal = numberRows_ + numberColumns;
220    CoinWorkDouble *region1Save = new CoinWorkDouble[numberTotal];
221    for (iColumn = 0; iColumn < numberTotal; iColumn++) {
222      region1[iColumn] *= diagonal[iColumn];
223      region1Save[iColumn] = region1[iColumn];
224    }
225    multiplyAdd(region1 + numberColumns, numberRows_, -1.0, region2, 1.0);
226    model_->clpMatrix()->times(1.0, region1, region2);
227    CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
228    CoinWorkDouble scale = 1.0;
229    CoinWorkDouble unscale = 1.0;
230    if (maximumRHS > 1.0e-30) {
231      if (maximumRHS <= 0.5) {
232        CoinWorkDouble factor = 2.0;
233        while (maximumRHS <= 0.5) {
234          maximumRHS *= factor;
235          scale *= factor;
236        } /* endwhile */
237      } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
238        CoinWorkDouble factor = 0.5;
239        while (maximumRHS >= 2.0) {
240          maximumRHS *= factor;
241          scale *= factor;
242        } /* endwhile */
243      }
244      unscale = diagonalScaleFactor / scale;
245    } else {
246      //effectively zero
247      scale = 0.0;
248      unscale = 0.0;
249    }
250    multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
251    solve(region2);
252    multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
253    multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns, 0.0);
254    CoinZeroN(region1, numberColumns);
255    model_->clpMatrix()->transposeTimes(1.0, region2, region1);
256    for (iColumn = 0; iColumn < numberTotal; iColumn++)
257      region1[iColumn] = region1[iColumn] * diagonal[iColumn] - region1Save[iColumn];
258    delete[] region1Save;
259  } else {
260    // KKT
261    int numberRowsModel = model_->numberRows();
262    int numberColumns = model_->numberColumns();
263    int numberTotal = numberColumns + numberRowsModel;
264    CoinWorkDouble *array = new CoinWorkDouble[numberRows_];
265    CoinMemcpyN(region1, numberTotal, array);
266    CoinMemcpyN(region2, numberRowsModel, array + numberTotal);
267    assert(numberRows_ >= numberRowsModel + numberTotal);
268    solve(array);
269    int iRow;
270    for (iRow = 0; iRow < numberTotal; iRow++) {
271      if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
272        COIN_DETAIL_PRINT(printf("row region1 %d dropped %g\n", iRow, array[iRow]));
273      }
274    }
275    for (; iRow < numberRows_; iRow++) {
276      if (rowsDropped_[iRow] && CoinAbs(array[iRow]) > 1.0e-8) {
277        COIN_DETAIL_PRINT(printf("row region2 %d dropped %g\n", iRow, array[iRow]));
278      }
279    }
280    CoinMemcpyN(array + numberTotal, numberRowsModel, region2);
281    CoinMemcpyN(array, numberTotal, region1);
282    delete[] array;
283  }
284}
285//-------------------------------------------------------------------
286// Clone
287//-------------------------------------------------------------------
288ClpCholeskyBase *ClpCholeskyBase::clone() const
289{
290  return new ClpCholeskyBase(*this);
291}
292// Forms ADAT - returns nonzero if not enough memory
293int ClpCholeskyBase::preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT)
294{
295  delete rowCopy_;
296  rowCopy_ = model_->clpMatrix()->reverseOrderedCopy();
297  if (!doKKT) {
298    numberRows_ = model_->numberRows();
299    rowsDropped_ = new char[numberRows_];
300    memset(rowsDropped_, 0, numberRows_);
301    numberRowsDropped_ = 0;
302    // Space for starts
303    choleskyStart_ = new int[numberRows_ + 1];
304    const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
305    const int *columnLength = model_->clpMatrix()->getVectorLengths();
306    const int *row = model_->clpMatrix()->getIndices();
307    const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
308    const int *rowLength = rowCopy_->getVectorLengths();
309    const int *column = rowCopy_->getIndices();
310    // We need two arrays for counts
311    int *which = new int[numberRows_];
312    int *used = new int[numberRows_ + 1];
313    CoinZeroN(used, numberRows_);
314    int iRow;
315    sizeFactor_ = 0;
316    int numberColumns = model_->numberColumns();
317    int numberDense = 0;
318    //denseThreshold_=3;
319    if (denseThreshold_ > 0) {
320      delete[] whichDense_;
321      delete[] denseColumn_;
322      delete dense_;
323      whichDense_ = new char[numberColumns];
324      int iColumn;
325      used[numberRows_] = 0;
326      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
327        int length = columnLength[iColumn];
328        used[length] += 1;
329      }
330      int nLong = 0;
331      int stop = CoinMax(denseThreshold_ / 2, 100);
332      for (iRow = numberRows_; iRow >= stop; iRow--) {
333        if (used[iRow])
334          COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
335        nLong += used[iRow];
336        if (nLong > 50 || nLong > (numberColumns >> 2))
337          break;
338      }
339      CoinZeroN(used, numberRows_);
340      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
341        if (columnLength[iColumn] < denseThreshold_) {
342          whichDense_[iColumn] = 0;
343        } else {
344          whichDense_[iColumn] = 1;
345          numberDense++;
346        }
347      }
348      if (!numberDense || numberDense > 100) {
349        // free
350        delete[] whichDense_;
351        whichDense_ = NULL;
352        denseColumn_ = NULL;
353        dense_ = NULL;
354      } else {
355        // space for dense columns
356        denseColumn_ = new longDouble[numberDense * numberRows_];
357        // dense cholesky
358        dense_ = new ClpCholeskyDense();
359        dense_->reserveSpace(NULL, numberDense);
360        COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
361      }
362    }
363    int offset = includeDiagonal ? 0 : 1;
364    if (lowerTriangular)
365      offset = -offset;
366    for (iRow = 0; iRow < numberRows_; iRow++) {
367      int number = 0;
368      // make sure diagonal exists if includeDiagonal
369      if (!offset) {
370        which[0] = iRow;
371        used[iRow] = 1;
372        number = 1;
373      }
374      CoinBigIndex startRow = rowStart[iRow];
375      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
376      if (lowerTriangular) {
377        for (CoinBigIndex k = startRow; k < endRow; k++) {
378          int iColumn = column[k];
379          if (!whichDense_ || !whichDense_[iColumn]) {
380            CoinBigIndex start = columnStart[iColumn];
381            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
382            for (CoinBigIndex j = start; j < end; j++) {
383              int jRow = row[j];
384              if (jRow <= iRow + offset) {
385                if (!used[jRow]) {
386                  used[jRow] = 1;
387                  which[number++] = jRow;
388                }
389              }
390            }
391          }
392        }
393      } else {
394        for (CoinBigIndex k = startRow; k < endRow; k++) {
395          int iColumn = column[k];
396          if (!whichDense_ || !whichDense_[iColumn]) {
397            CoinBigIndex start = columnStart[iColumn];
398            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
399            for (CoinBigIndex j = start; j < end; j++) {
400              int jRow = row[j];
401              if (jRow >= iRow + offset) {
402                if (!used[jRow]) {
403                  used[jRow] = 1;
404                  which[number++] = jRow;
405                }
406              }
407            }
408          }
409        }
410      }
411      sizeFactor_ += number;
412      int j;
413      for (j = 0; j < number; j++)
414        used[which[j]] = 0;
415    }
416    delete[] which;
417    // Now we have size - create arrays and fill in
418    try {
419      choleskyRow_ = new int[sizeFactor_];
420    } catch (...) {
421      // no memory
422      delete[] choleskyStart_;
423      choleskyStart_ = NULL;
424      return -1;
425    }
426    sizeFactor_ = 0;
427    which = choleskyRow_;
428    for (iRow = 0; iRow < numberRows_; iRow++) {
429      int number = 0;
430      // make sure diagonal exists if includeDiagonal
431      if (!offset) {
432        which[0] = iRow;
433        used[iRow] = 1;
434        number = 1;
435      }
436      choleskyStart_[iRow] = sizeFactor_;
437      CoinBigIndex startRow = rowStart[iRow];
438      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
439      if (lowerTriangular) {
440        for (CoinBigIndex k = startRow; k < endRow; k++) {
441          int iColumn = column[k];
442          if (!whichDense_ || !whichDense_[iColumn]) {
443            CoinBigIndex start = columnStart[iColumn];
444            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
445            for (CoinBigIndex j = start; j < end; j++) {
446              int jRow = row[j];
447              if (jRow <= iRow + offset) {
448                if (!used[jRow]) {
449                  used[jRow] = 1;
450                  which[number++] = jRow;
451                }
452              }
453            }
454          }
455        }
456      } else {
457        for (CoinBigIndex k = startRow; k < endRow; k++) {
458          int iColumn = column[k];
459          if (!whichDense_ || !whichDense_[iColumn]) {
460            CoinBigIndex start = columnStart[iColumn];
461            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
462            for (CoinBigIndex j = start; j < end; j++) {
463              int jRow = row[j];
464              if (jRow >= iRow + offset) {
465                if (!used[jRow]) {
466                  used[jRow] = 1;
467                  which[number++] = jRow;
468                }
469              }
470            }
471          }
472        }
473      }
474      sizeFactor_ += number;
475      int j;
476      for (j = 0; j < number; j++)
477        used[which[j]] = 0;
478      // Sort
479      std::sort(which, which + number);
480      // move which on
481      which += number;
482    }
483    choleskyStart_[numberRows_] = sizeFactor_;
484    delete[] used;
485    return 0;
486  } else {
487    int numberRowsModel = model_->numberRows();
488    int numberColumns = model_->numberColumns();
489    int numberTotal = numberColumns + numberRowsModel;
490    numberRows_ = 2 * numberRowsModel + numberColumns;
491    rowsDropped_ = new char[numberRows_];
492    memset(rowsDropped_, 0, numberRows_);
493    numberRowsDropped_ = 0;
494    CoinPackedMatrix *quadratic = NULL;
495    ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
496    if (quadraticObj)
497      quadratic = quadraticObj->quadraticObjective();
498    CoinBigIndex numberElements = model_->clpMatrix()->getNumElements();
499    numberElements = numberElements + 2 * numberRowsModel + numberTotal;
500    if (quadratic)
501      numberElements += quadratic->getNumElements();
502    // Space for starts
503    choleskyStart_ = new int[numberRows_ + 1];
504    const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
505    const int *columnLength = model_->clpMatrix()->getVectorLengths();
506    const int *row = model_->clpMatrix()->getIndices();
507    //const double * element = model_->clpMatrix()->getElements();
508    // Now we have size - create arrays and fill in
509    try {
510      choleskyRow_ = new int[numberElements];
511    } catch (...) {
512      // no memory
513      delete[] choleskyStart_;
514      choleskyStart_ = NULL;
515      return -1;
516    }
517    int iRow, iColumn;
518
519    sizeFactor_ = 0;
520    // matrix
521    if (lowerTriangular) {
522      if (!quadratic) {
523        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
524          choleskyStart_[iColumn] = sizeFactor_;
525          choleskyRow_[sizeFactor_++] = iColumn;
526          CoinBigIndex start = columnStart[iColumn];
527          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
528          if (!includeDiagonal)
529            start++;
530          for (CoinBigIndex j = start; j < end; j++) {
531            choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
532          }
533        }
534      } else {
535        // Quadratic
536        const int *columnQuadratic = quadratic->getIndices();
537        const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
538        const int *columnQuadraticLength = quadratic->getVectorLengths();
539        //const double * quadraticElement = quadratic->getElements();
540        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
541          choleskyStart_[iColumn] = sizeFactor_;
542          if (includeDiagonal)
543            choleskyRow_[sizeFactor_++] = iColumn;
544          CoinBigIndex j;
545          for (j = columnQuadraticStart[iColumn];
546               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
547            int jColumn = columnQuadratic[j];
548            if (jColumn > iColumn)
549              choleskyRow_[sizeFactor_++] = jColumn;
550          }
551          CoinBigIndex start = columnStart[iColumn];
552          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
553          for (j = start; j < end; j++) {
554            choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
555          }
556        }
557      }
558      // slacks
559      for (; iColumn < numberTotal; iColumn++) {
560        choleskyStart_[iColumn] = sizeFactor_;
561        if (includeDiagonal)
562          choleskyRow_[sizeFactor_++] = iColumn;
563        choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
564      }
565      // Transpose - nonzero diagonal (may regularize)
566      for (iRow = 0; iRow < numberRowsModel; iRow++) {
567        choleskyStart_[iRow + numberTotal] = sizeFactor_;
568        // diagonal
569        if (includeDiagonal)
570          choleskyRow_[sizeFactor_++] = iRow + numberTotal;
571      }
572      choleskyStart_[numberRows_] = sizeFactor_;
573    } else {
574      // transpose
575      ClpMatrixBase *rowCopy = model_->clpMatrix()->reverseOrderedCopy();
576      const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
577      const int *rowLength = rowCopy->getVectorLengths();
578      const int *column = rowCopy->getIndices();
579      if (!quadratic) {
580        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
581          choleskyStart_[iColumn] = sizeFactor_;
582          if (includeDiagonal)
583            choleskyRow_[sizeFactor_++] = iColumn;
584        }
585      } else {
586        // Quadratic
587        // transpose
588        CoinPackedMatrix quadraticT;
589        quadraticT.reverseOrderedCopyOf(*quadratic);
590        const int *columnQuadratic = quadraticT.getIndices();
591        const CoinBigIndex *columnQuadraticStart = quadraticT.getVectorStarts();
592        const int *columnQuadraticLength = quadraticT.getVectorLengths();
593        //const double * quadraticElement = quadraticT.getElements();
594        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
595          choleskyStart_[iColumn] = sizeFactor_;
596          for (CoinBigIndex j = columnQuadraticStart[iColumn];
597               j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
598            int jColumn = columnQuadratic[j];
599            if (jColumn < iColumn)
600              choleskyRow_[sizeFactor_++] = jColumn;
601          }
602          if (includeDiagonal)
603            choleskyRow_[sizeFactor_++] = iColumn;
604        }
605      }
606      int iRow;
607      // slacks
608      for (iRow = 0; iRow < numberRowsModel; iRow++) {
609        choleskyStart_[iRow + numberColumns] = sizeFactor_;
610        if (includeDiagonal)
611          choleskyRow_[sizeFactor_++] = iRow + numberColumns;
612      }
613      for (iRow = 0; iRow < numberRowsModel; iRow++) {
614        choleskyStart_[iRow + numberTotal] = sizeFactor_;
615        CoinBigIndex start = rowStart[iRow];
616        CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
617        for (CoinBigIndex j = start; j < end; j++) {
618          choleskyRow_[sizeFactor_++] = column[j];
619        }
620        // slack
621        choleskyRow_[sizeFactor_++] = numberColumns + iRow;
622        if (includeDiagonal)
623          choleskyRow_[sizeFactor_++] = iRow + numberTotal;
624      }
625      choleskyStart_[numberRows_] = sizeFactor_;
626    }
627  }
628  return 0;
629}
630/* Orders rows and saves pointer to matrix.and model */
631int ClpCholeskyBase::order(ClpInterior *model)
632{
633  model_ = model;
634#define BASE_ORDER 2
635#if BASE_ORDER > 0
636  if (model_->numberRows() > 6) {
637    if (preOrder(doKKT_, true, doKKT_))
638      return -1;
639    //rowsDropped_ = new char [numberRows_];
640    numberRowsDropped_ = 0;
641    memset(rowsDropped_, 0, numberRows_);
642    //rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
643    // approximate minimum degree
644    return orderAMD();
645  }
646#endif
647  int numberRowsModel = model_->numberRows();
648  int numberColumns = model_->numberColumns();
649  int numberTotal = numberColumns + numberRowsModel;
650  CoinPackedMatrix *quadratic = NULL;
651  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
652  if (quadraticObj)
653    quadratic = quadraticObj->quadraticObjective();
654  if (!doKKT_) {
655    numberRows_ = model->numberRows();
656  } else {
657    numberRows_ = 2 * numberRowsModel + numberColumns;
658  }
659  rowsDropped_ = new char[numberRows_];
660  numberRowsDropped_ = 0;
661  memset(rowsDropped_, 0, numberRows_);
662  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
663  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
664  const int *columnLength = model_->clpMatrix()->getVectorLengths();
665  const int *row = model_->clpMatrix()->getIndices();
666  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
667  const int *rowLength = rowCopy_->getVectorLengths();
668  const int *column = rowCopy_->getIndices();
669  // We need arrays for counts
670  int *which = new int[numberRows_];
671  int *used = new int[numberRows_ + 1];
672  int *count = new int[numberRows_];
673  CoinZeroN(count, numberRows_);
674  CoinZeroN(used, numberRows_);
675  int iRow;
676  sizeFactor_ = 0;
677  permute_ = new int[numberRows_];
678  for (iRow = 0; iRow < numberRows_; iRow++)
679    permute_[iRow] = iRow;
680  if (!doKKT_) {
681    int numberDense = 0;
682    if (denseThreshold_ > 0) {
683      delete[] whichDense_;
684      delete[] denseColumn_;
685      delete dense_;
686      whichDense_ = new char[numberColumns];
687      int iColumn;
688      used[numberRows_] = 0;
689      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
690        int length = columnLength[iColumn];
691        used[length] += 1;
692      }
693      int nLong = 0;
694      int stop = CoinMax(denseThreshold_ / 2, 100);
695      for (iRow = numberRows_; iRow >= stop; iRow--) {
696        if (used[iRow])
697          COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow));
698        nLong += used[iRow];
699        if (nLong > 50 || nLong > (numberColumns >> 2))
700          break;
701      }
702      CoinZeroN(used, numberRows_);
703      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
704        if (columnLength[iColumn] < denseThreshold_) {
705          whichDense_[iColumn] = 0;
706        } else {
707          whichDense_[iColumn] = 1;
708          numberDense++;
709        }
710      }
711      if (!numberDense || numberDense > 100) {
712        // free
713        delete[] whichDense_;
714        whichDense_ = NULL;
715        denseColumn_ = NULL;
716        dense_ = NULL;
717      } else {
718        // space for dense columns
719        denseColumn_ = new longDouble[numberDense * numberRows_];
720        // dense cholesky
721        dense_ = new ClpCholeskyDense();
722        dense_->reserveSpace(NULL, numberDense);
723        COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense));
724      }
725    }
726    /*
727             Get row counts and size
728          */
729    for (iRow = 0; iRow < numberRows_; iRow++) {
730      int number = 1;
731      // make sure diagonal exists
732      which[0] = iRow;
733      used[iRow] = 1;
734      CoinBigIndex startRow = rowStart[iRow];
735      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
736      for (CoinBigIndex k = startRow; k < endRow; k++) {
737        int iColumn = column[k];
738        if (!whichDense_ || !whichDense_[iColumn]) {
739          CoinBigIndex start = columnStart[iColumn];
740          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
741          for (CoinBigIndex j = start; j < end; j++) {
742            int jRow = row[j];
743            if (jRow < iRow) {
744              if (!used[jRow]) {
745                used[jRow] = 1;
746                which[number++] = jRow;
747                count[jRow]++;
748              }
749            }
750          }
751        }
752      }
753      sizeFactor_ += number;
754      count[iRow] += number;
755      int j;
756      for (j = 0; j < number; j++)
757        used[which[j]] = 0;
758    }
759    CoinSort_2(count, count + numberRows_, permute_);
760  } else {
761    // KKT
762    CoinBigIndex numberElements = model_->clpMatrix()->getNumElements();
763    numberElements = numberElements + 2 * numberRowsModel + numberTotal;
764    if (quadratic)
765      numberElements += quadratic->getNumElements();
766    // off diagonal
767    numberElements -= numberRows_;
768    sizeFactor_ = static_cast< int >(numberElements);
769    // If we sort we need to redo symbolic etc
770  }
771  delete[] which;
772  delete[] used;
773  delete[] count;
774  permuteInverse_ = new int[numberRows_];
775  for (iRow = 0; iRow < numberRows_; iRow++) {
776    //permute_[iRow]=iRow; // force no permute
777    //permute_[iRow]=numberRows_-1-iRow; // force odd permute
778    //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
779    permuteInverse_[permute_[iRow]] = iRow;
780  }
781  return 0;
782}
783#if BASE_ORDER == 1
784/* Orders rows and saves pointer to matrix.and model */
785int ClpCholeskyBase::orderAMD()
786{
787  permuteInverse_ = new int[numberRows_];
788  permute_ = new int[numberRows_];
789  // Do ordering
790  int returnCode = 0;
791  // get more space and full matrix
792  int space = 6 * sizeFactor_ + 100000;
793  int *temp = new int[space];
794  int *which = new int[2 * numberRows_];
795  CoinBigIndex *tempStart = new CoinBigIndex[numberRows_ + 1];
796  memset(which, 0, numberRows_ * sizeof(int));
797  for (int iRow = 0; iRow < numberRows_; iRow++) {
798    which[iRow] += choleskyStart_[iRow + 1] - choleskyStart_[iRow] - 1;
799    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
800      int jRow = choleskyRow_[j];
801      which[jRow]++;
802    }
803  }
804  CoinBigIndex sizeFactor = 0;
805  for (int iRow = 0; iRow < numberRows_; iRow++) {
806    int length = which[iRow];
807    permute_[iRow] = length;
808    tempStart[iRow] = sizeFactor;
809    which[iRow] = sizeFactor;
810    sizeFactor += length;
811  }
812  for (int iRow = 0; iRow < numberRows_; iRow++) {
813    assert(choleskyRow_[choleskyStart_[iRow]] == iRow);
814    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
815      int jRow = choleskyRow_[j];
816      int put = which[iRow];
817      temp[put] = jRow;
818      which[iRow]++;
819      put = which[jRow];
820      temp[put] = iRow;
821      which[jRow]++;
822    }
823  }
824  for (int iRow = 1; iRow < numberRows_; iRow++)
825    assert(which[iRow - 1] == tempStart[iRow]);
826  CoinBigIndex lastSpace = sizeFactor;
827  delete[] choleskyRow_;
828  choleskyRow_ = temp;
829  delete[] choleskyStart_;
830  choleskyStart_ = tempStart;
831  // linked lists of sizes and lengths
832  int *first = new int[numberRows_];
833  int *next = new int[numberRows_];
834  int *previous = new int[numberRows_];
835  char *mark = new char[numberRows_];
836  memset(mark, 0, numberRows_);
837  CoinBigIndex *sort = new CoinBigIndex[numberRows_];
838  int *order = new int[numberRows_];
839  for (int iRow = 0; iRow < numberRows_; iRow++) {
840    first[iRow] = -1;
841    next[iRow] = -1;
842    previous[iRow] = -1;
843    permuteInverse_[iRow] = -1;
844  }
845  int large = 1000 + 2 * numberRows_;
846  for (int iRow = 0; iRow < numberRows_; iRow++) {
847    int n = permute_[iRow];
848    if (first[n] < 0) {
849      first[n] = iRow;
850      previous[iRow] = n + large;
851      next[iRow] = n + 2 * large;
852    } else {
853      int k = first[n];
854      assert(k < numberRows_);
855      first[n] = iRow;
856      previous[iRow] = n + large;
857      assert(previous[k] == n + large);
858      next[iRow] = k;
859      previous[k] = iRow;
860    }
861  }
862  int smallest = 0;
863  int done = 0;
864  int numberCompressions = 0;
865  int numberExpansions = 0;
866  while (smallest < numberRows_) {
867    if (first[smallest] < 0 || first[smallest] > numberRows_) {
868      smallest++;
869      continue;
870    }
871    int iRow = first[smallest];
872    if (false) {
873      int kRow = -1;
874      int ss = 999999;
875      for (int jRow = numberRows_ - 1; jRow >= 0; jRow--) {
876        if (permuteInverse_[jRow] < 0) {
877          int length = permute_[jRow];
878          assert(length > 0);
879          for (CoinBigIndex j = choleskyStart_[jRow];
880               j < choleskyStart_[jRow] + length; j++) {
881            int jjRow = choleskyRow_[j];
882            assert(permuteInverse_[jjRow] < 0);
883          }
884          if (length < ss) {
885            ss = length;
886            kRow = jRow;
887          }
888        }
889      }
890      assert(smallest == ss);
891      printf("list chose %d - full chose %d - length %d\n",
892        iRow, kRow, ss);
893    }
894    int kNext = next[iRow];
895    first[smallest] = kNext;
896    if (kNext < numberRows_)
897      previous[kNext] = previous[iRow];
898    previous[iRow] = -1;
899    next[iRow] = -1;
900    permuteInverse_[iRow] = done;
901    done++;
902    // Now add edges
903    CoinBigIndex start = choleskyStart_[iRow];
904    CoinBigIndex end = choleskyStart_[iRow] + permute_[iRow];
905    int nSave = 0;
906    for (CoinBigIndex k = start; k < end; k++) {
907      int kRow = choleskyRow_[k];
908      assert(permuteInverse_[kRow] < 0);
909      //if (permuteInverse_[kRow]<0)
910      which[nSave++] = kRow;
911    }
912    for (int i = 0; i < nSave; i++) {
913      int kRow = which[i];
914      int length = permute_[kRow];
915      CoinBigIndex start = choleskyStart_[kRow];
916      CoinBigIndex end = choleskyStart_[kRow] + length;
917      for (CoinBigIndex j = start; j < end; j++) {
918        int jRow = choleskyRow_[j];
919        mark[jRow] = 1;
920      }
921      mark[kRow] = 1;
922      int n = nSave;
923      for (int j = 0; j < nSave; j++) {
924        int jRow = which[j];
925        if (!mark[jRow]) {
926          which[n++] = jRow;
927        }
928      }
929      for (CoinBigIndex j = start; j < end; j++) {
930        int jRow = choleskyRow_[j];
931        mark[jRow] = 0;
932      }
933      mark[kRow] = 0;
934      if (n > nSave) {
935        bool inPlace = (n - nSave == 1);
936        if (!inPlace) {
937          // extend
938          int length = n - nSave + end - start;
939          if (length + lastSpace > space) {
940            // need to compress
941            numberCompressions++;
942            int newN = 0;
943            for (int iRow = 0; iRow < numberRows_; iRow++) {
944              CoinBigIndex start = choleskyStart_[iRow];
945              if (permuteInverse_[iRow] < 0) {
946                sort[newN] = start;
947                order[newN++] = iRow;
948              } else {
949                choleskyStart_[iRow] = -1;
950                permute_[iRow] = 0;
951              }
952            }
953            CoinSort_2(sort, sort + newN, order);
954            sizeFactor = 0;
955            for (int k = 0; k < newN; k++) {
956              int iRow = order[k];
957              int length = permute_[iRow];
958              CoinBigIndex start = choleskyStart_[iRow];
959              choleskyStart_[iRow] = sizeFactor;
960              for (int j = 0; j < length; j++)
961                choleskyRow_[sizeFactor + j] = choleskyRow_[start + j];
962              sizeFactor += length;
963            }
964            lastSpace = sizeFactor;
965            if ((sizeFactor + numberRows_ + 1000) * 4 > 3 * space) {
966              // need to expand
967              numberExpansions++;
968              space = (3 * space) / 2;
969              int *temp = new int[space];
970              memcpy(temp, choleskyRow_, sizeFactor * sizeof(int));
971              delete[] choleskyRow_;
972              choleskyRow_ = temp;
973            }
974          }
975        }
976        // Now add
977        start = choleskyStart_[kRow];
978        end = choleskyStart_[kRow] + permute_[kRow];
979        if (!inPlace)
980          choleskyStart_[kRow] = lastSpace;
981        CoinBigIndex put = choleskyStart_[kRow];
982        for (CoinBigIndex j = start; j < end; j++) {
983          int jRow = choleskyRow_[j];
984          if (permuteInverse_[jRow] < 0)
985            choleskyRow_[put++] = jRow;
986        }
987        for (int j = nSave; j < n; j++) {
988          int jRow = which[j];
989          choleskyRow_[put++] = jRow;
990        }
991        if (!inPlace) {
992          permute_[kRow] = put - lastSpace;
993          lastSpace = put;
994        } else {
995          permute_[kRow] = put - choleskyStart_[kRow];
996        }
997      } else {
998        // pack down for new counts
999        CoinBigIndex put = start;
1000        for (CoinBigIndex j = start; j < end; j++) {
1001          int jRow = choleskyRow_[j];
1002          if (permuteInverse_[jRow] < 0)
1003            choleskyRow_[put++] = jRow;
1004        }
1005        permute_[kRow] = put - start;
1006      }
1007      // take out
1008      int iNext = next[kRow];
1009      int iPrevious = previous[kRow];
1010      if (iPrevious < numberRows_) {
1011        next[iPrevious] = iNext;
1012      } else {
1013        assert(iPrevious == length + large);
1014        first[length] = iNext;
1015      }
1016      if (iNext < numberRows_) {
1017        previous[iNext] = iPrevious;
1018      } else {
1019        assert(iNext == length + 2 * large);
1020      }
1021      // put in
1022      length = permute_[kRow];
1023      smallest = CoinMin(smallest, length);
1024      if (first[length] < 0 || first[length] > numberRows_) {
1025        first[length] = kRow;
1026        previous[kRow] = length + large;
1027        next[kRow] = length + 2 * large;
1028      } else {
1029        int k = first[length];
1030        assert(k < numberRows_);
1031        first[length] = kRow;
1032        previous[kRow] = length + large;
1033        assert(previous[k] == length + large);
1034        next[kRow] = k;
1035        previous[k] = kRow;
1036      }
1037    }
1038  }
1039  // do rest
1040  for (int iRow = 0; iRow < numberRows_; iRow++) {
1041    if (permuteInverse_[iRow] < 0)
1042      permuteInverse_[iRow] = done++;
1043  }
1044  COIN_DETAIL_PRINT(printf("%d compressions, %d expansions\n",
1045    numberCompressions, numberExpansions));
1046  assert(done == numberRows_);
1047  delete[] sort;
1048  delete[] order;
1049  delete[] which;
1050  delete[] mark;
1051  delete[] first;
1052  delete[] next;
1053  delete[] previous;
1054  delete[] choleskyRow_;
1055  choleskyRow_ = NULL;
1056  delete[] choleskyStart_;
1057  choleskyStart_ = NULL;
1058#ifndef NDEBUG
1059  for (int iRow = 0; iRow < numberRows_; iRow++) {
1060    permute_[iRow] = -1;
1061  }
1062#endif
1063  for (int iRow = 0; iRow < numberRows_; iRow++) {
1064    permute_[permuteInverse_[iRow]] = iRow;
1065  }
1066#ifndef NDEBUG
1067  for (int iRow = 0; iRow < numberRows_; iRow++) {
1068    assert(permute_[iRow] >= 0 && permute_[iRow] < numberRows_);
1069  }
1070#endif
1071  return returnCode;
1072}
1073#elif BASE_ORDER == 2
1074/*----------------------------------------------------------------------------*/
1075/*      (C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
1076/*----------------------------------------------------------------------------*/
1077
1078/*  A compact no-frills Approximate Minimum Local Fill ordering code.
1079
1080    References:
1081
1082[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
1083    Edward Rothberg, SGI Manuscript, April 1996.
1084[2] An Approximate Minimum Degree Ordering Algorithm.
1085    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
1086    University of Florida, December 1994.
1087*/
1088
1089#include <math.h>
1090#include <stdlib.h>
1091
1092typedef int WSI;
1093
1094#define NORDTHRESH 7
1095#define MAXIW 2147000000
1096
1097//#define WSSMP_DBG
1098#ifdef WSSMP_DBG
1099static void chk(WSI ind, WSI i, WSI lim)
1100{
1101
1102  if (i <= 0 || i > lim) {
1103    printf("%d: chk: myamlf: WAR**: i, lim = %d, %d\n", ind, i, lim);
1104  }
1105}
1106#endif
1107
1108static void
1109myamlf(WSI n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI varbl[],
1110  WSI snxt[], WSI perm[], WSI invp[], WSI head[], WSI lsize[],
1111  WSI flag[], WSI erscore[], WSI locaux, WSI adjln, WSI speed)
1112{
1113  WSI l, i, j, k;
1114  double x, y;
1115  WSI maxmum, fltag, nodeg, scln, nm1, deg, tn,
1116    locatns, ipp, jpp, nnode, numpiv, node,
1117    nodeln, currloc, counter, numii, mindeg,
1118    i0, i1, i2, i4, i5, i6, i7, i9,
1119    j0, j1, j2, j3, j4, j5, j6, j7, j8, j9;
1120  /*                                             n cannot be less than NORDTHRESH
1121      if (n < 3) {
1122         if (n > 0) perm[0] = invp[0] = 1;
1123         if (n > 1) perm[1] = invp[1] = 2;
1124         return;
1125      }
1126     */
1127#ifdef WSSMP_DBG
1128  printf("Myamlf: n, locaux, adjln, speed = %d,%d,%d,%d\n", n, locaux, adjln, speed);
1129  for (i = 0; i < n; i++)
1130    flag[i] = 0;
1131  k = xadj[0];
1132  for (i = 1; i <= n; i++) {
1133    j = k + dgree[i - 1];
1134    if (j < k || k < 1)
1135      printf("ERR**: myamlf: i, j, k = %d,%d,%d\n", i, j, k);
1136    for (l = k; l < j; l++) {
1137      if (adjncy[l - 1] < 1 || adjncy[l - 1] > n || adjncy[l - 1] == i)
1138        printf("ERR**: myamlf: i, l, adjj[l] = %d,%d,%d\n", i, l, adjncy[l - 1]);
1139      if (flag[adjncy[l - 1] - 1] == i)
1140        printf("ERR**: myamlf: duplicate entry: %d - %d\n", i, adjncy[l - 1]);
1141      flag[adjncy[l - 1] - 1] = i;
1142      nm1 = adjncy[l - 1] - 1;
1143      for (fltag = xadj[nm1]; fltag < xadj[nm1] + dgree[nm1]; fltag++) {
1144        if (adjncy[fltag - 1] == i)
1145          goto L100;
1146      }
1147      printf("ERR**: Unsymmetric %d %d\n", i, fltag);
1148    L100:;
1149    }
1150    k = j;
1151    if (k > locaux)
1152      printf("ERR**: i, j, k, locaux = %d, %d, %d, %d\n",
1153        i, j, k, locaux);
1154  }
1155  if (n * (n - 1) < locaux - 1)
1156    printf("ERR**: myamlf: too many nnz, n, ne = %d, %d\n",
1157      n, adjln);
1158  for (i = 0; i < n; i++)
1159    flag[i] = 1;
1160  if (n > 10000)
1161    printf("Finished error checking in AMF\n");
1162  for (i = locaux; i <= adjln; i++)
1163    adjncy[i - 1] = -22;
1164#endif
1165
1166  maxmum = 0;
1167  mindeg = 1;
1168  fltag = 2;
1169  locatns = locaux - 1;
1170  nm1 = n - 1;
1171  counter = 1;
1172  for (l = 0; l < n; l++) {
1173    j = erscore[l];
1174#ifdef WSSMP_DBG
1175    chk(1, j, n);
1176#endif
1177    if (j > 0) {
1178      nnode = head[j - 1];
1179      if (nnode)
1180        perm[nnode - 1] = l + 1;
1181      snxt[l] = nnode;
1182      head[j - 1] = l + 1;
1183    } else {
1184      invp[l] = -(counter++);
1185      flag[l] = xadj[l] = 0;
1186    }
1187  }
1188  while (counter <= n) {
1189    for (deg = mindeg; head[deg - 1] < 1; deg++) {
1190    };
1191    nodeg = 0;
1192#ifdef WSSMP_DBG
1193    chk(2, deg, n - 1);
1194#endif
1195    node = head[deg - 1];
1196#ifdef WSSMP_DBG
1197    chk(3, node, n);
1198#endif
1199    mindeg = deg;
1200    nnode = snxt[node - 1];
1201    if (nnode)
1202      perm[nnode - 1] = 0;
1203    head[deg - 1] = nnode;
1204    nodeln = invp[node - 1];
1205    numpiv = varbl[node - 1];
1206    invp[node - 1] = -counter;
1207    counter += numpiv;
1208    varbl[node - 1] = -numpiv;
1209    if (nodeln) {
1210      j4 = locaux;
1211      i5 = lsize[node - 1] - nodeln;
1212      i2 = nodeln + 1;
1213      l = xadj[node - 1];
1214      for (i6 = 1; i6 <= i2; ++i6) {
1215        if (i6 == i2) {
1216          tn = node, i0 = l, scln = i5;
1217        } else {
1218#ifdef WSSMP_DBG
1219          chk(4, l, adjln);
1220#endif
1221          tn = adjncy[l - 1];
1222          l++;
1223#ifdef WSSMP_DBG
1224          chk(5, tn, n);
1225#endif
1226          i0 = xadj[tn - 1];
1227          scln = lsize[tn - 1];
1228        }
1229        for (i7 = 1; i7 <= scln; ++i7) {
1230#ifdef WSSMP_DBG
1231          chk(6, i0, adjln);
1232#endif
1233          i = adjncy[i0 - 1];
1234          i0++;
1235#ifdef WSSMP_DBG
1236          chk(7, i, n);
1237#endif
1238          numii = varbl[i - 1];
1239          if (numii > 0) {
1240            if (locaux > adjln) {
1241              lsize[node - 1] -= i6;
1242              xadj[node - 1] = l;
1243              if (!lsize[node - 1])
1244                xadj[node - 1] = 0;
1245              xadj[tn - 1] = i0;
1246              lsize[tn - 1] = scln - i7;
1247              if (!lsize[tn - 1])
1248                xadj[tn - 1] = 0;
1249              for (j = 1; j <= n; j++) {
1250                i4 = xadj[j - 1];
1251                if (i4 > 0) {
1252                  xadj[j - 1] = adjncy[i4 - 1];
1253#ifdef WSSMP_DBG
1254                  chk(8, i4, adjln);
1255#endif
1256                  adjncy[i4 - 1] = -j;
1257                }
1258              }
1259              i9 = j4 - 1;
1260              j6 = j7 = 1;
1261              while (j6 <= i9) {
1262#ifdef WSSMP_DBG
1263                chk(9, j6, adjln);
1264#endif
1265                j = -adjncy[j6 - 1];
1266                j6++;
1267                if (j > 0) {
1268#ifdef WSSMP_DBG
1269                  chk(10, j7, adjln);
1270                  chk(11, j, n);
1271#endif
1272                  adjncy[j7 - 1] = xadj[j - 1];
1273                  xadj[j - 1] = j7++;
1274                  j8 = lsize[j - 1] - 1 + j7;
1275                  for (; j7 < j8; j7++, j6++)
1276                    adjncy[j7 - 1] = adjncy[j6 - 1];
1277                }
1278              }
1279              for (j0 = j7; j4 < locaux; j4++, j7++) {
1280#ifdef WSSMP_DBG
1281                chk(12, j4, adjln);
1282#endif
1283                adjncy[j7 - 1] = adjncy[j4 - 1];
1284              }
1285              j4 = j0;
1286              locaux = j7;
1287              i0 = xadj[tn - 1];
1288              l = xadj[node - 1];
1289            }
1290            adjncy[locaux - 1] = i;
1291            locaux++;
1292            varbl[i - 1] = -numii;
1293            nodeg += numii;
1294            ipp = perm[i - 1];
1295            nnode = snxt[i - 1];
1296#ifdef WSSMP_DBG
1297            if (ipp)
1298              chk(13, ipp, n);
1299            if (nnode)
1300              chk(14, nnode, n);
1301#endif
1302            if (ipp)
1303              snxt[ipp - 1] = nnode;
1304            else
1305              head[erscore[i - 1] - 1] = nnode;
1306            if (nnode)
1307              perm[nnode - 1] = ipp;
1308          }
1309        }
1310        if (tn != node) {
1311          flag[tn - 1] = 0, xadj[tn - 1] = -node;
1312        }
1313      }
1314      currloc = (j5 = locaux) - j4;
1315      locatns += currloc;
1316    } else {
1317      i1 = (j4 = xadj[node - 1]) + lsize[node - 1];
1318      for (j = j5 = j4; j < i1; j++) {
1319#ifdef WSSMP_DBG
1320        chk(15, j, adjln);
1321#endif
1322        i = adjncy[j - 1];
1323#ifdef WSSMP_DBG
1324        chk(16, i, n);
1325#endif
1326        numii = varbl[i - 1];
1327        if (numii > 0) {
1328          nodeg += numii;
1329          varbl[i - 1] = -numii;
1330#ifdef WSSMP_DBG
1331          chk(17, j5, adjln);
1332#endif
1333          adjncy[j5 - 1] = i;
1334          ipp = perm[i - 1];
1335          nnode = snxt[i - 1];
1336          j5++;
1337#ifdef WSSMP_DBG
1338          if (ipp)
1339            chk(18, ipp, n);
1340          if (nnode)
1341            chk(19, nnode, n);
1342#endif
1343          if (ipp)
1344            snxt[ipp - 1] = nnode;
1345          else
1346            head[erscore[i - 1] - 1] = nnode;
1347          if (nnode)
1348            perm[nnode - 1] = ipp;
1349        }
1350      }
1351      currloc = 0;
1352    }
1353#ifdef WSSMP_DBG
1354    chk(20, node, n);
1355#endif
1356    lsize[node - 1] = j5 - (xadj[node - 1] = j4);
1357    dgree[node - 1] = nodeg;
1358    if (fltag + n < 0 || fltag + n > MAXIW) {
1359      for (i = 1; i <= n; i++)
1360        if (flag[i - 1])
1361          flag[i - 1] = 1;
1362      fltag = 2;
1363    }
1364    for (j3 = j4; j3 < j5; j3++) {
1365#ifdef WSSMP_DBG
1366      chk(21, j3, adjln);
1367#endif
1368      i = adjncy[j3 - 1];
1369#ifdef WSSMP_DBG
1370      chk(22, i, n);
1371#endif
1372      j = invp[i - 1];
1373      if (j > 0) {
1374        i4 = fltag - (numii = -varbl[i - 1]);
1375        i2 = xadj[i - 1] + j;
1376        for (l = xadj[i - 1]; l < i2; l++) {
1377#ifdef WSSMP_DBG
1378          chk(23, l, adjln);
1379#endif
1380          tn = adjncy[l - 1];
1381#ifdef WSSMP_DBG
1382          chk(24, tn, n);
1383#endif
1384          j9 = flag[tn - 1];
1385          if (j9 >= fltag)
1386            j9 -= numii;
1387          else if (j9)
1388            j9 = dgree[tn - 1] + i4;
1389          flag[tn - 1] = j9;
1390        }
1391      }
1392    }
1393    for (j3 = j4; j3 < j5; j3++) {
1394#ifdef WSSMP_DBG
1395      chk(25, j3, adjln);
1396#endif
1397      i = adjncy[j3 - 1];
1398      i5 = deg = 0;
1399#ifdef WSSMP_DBG
1400      chk(26, i, n);
1401#endif
1402      j1 = (i4 = xadj[i - 1]) + invp[i - 1];
1403      for (l = j0 = i4; l < j1; l++) {
1404#ifdef WSSMP_DBG
1405        chk(27, l, adjln);
1406#endif
1407        tn = adjncy[l - 1];
1408#ifdef WSSMP_DBG
1409        chk(70, tn, n);
1410#endif
1411        j8 = flag[tn - 1];
1412        if (j8) {
1413          deg += (j8 - fltag);
1414#ifdef WSSMP_DBG
1415          chk(28, i4, adjln);
1416#endif
1417          adjncy[i4 - 1] = tn;
1418          i5 += tn;
1419          i4++;
1420          while (i5 >= nm1)
1421            i5 -= nm1;
1422        }
1423      }
1424      invp[i - 1] = (j2 = i4) - j0 + 1;
1425      i2 = j0 + lsize[i - 1];
1426      for (l = j1; l < i2; l++) {
1427#ifdef WSSMP_DBG
1428        chk(29, l, adjln);
1429#endif
1430        j = adjncy[l - 1];
1431#ifdef WSSMP_DBG
1432        chk(30, j, n);
1433#endif
1434        numii = varbl[j - 1];
1435        if (numii > 0) {
1436          deg += numii;
1437#ifdef WSSMP_DBG
1438          chk(31, i4, adjln);
1439#endif
1440          adjncy[i4 - 1] = j;
1441          i5 += j;
1442          i4++;
1443          while (i5 >= nm1)
1444            i5 -= nm1;
1445        }
1446      }
1447      if (invp[i - 1] == 1 && j2 == i4) {
1448        numii = -varbl[i - 1];
1449        xadj[i - 1] = -node;
1450        nodeg -= numii;
1451        counter += numii;
1452        numpiv += numii;
1453        invp[i - 1] = varbl[i - 1] = 0;
1454      } else {
1455        if (dgree[i - 1] > deg)
1456          dgree[i - 1] = deg;
1457#ifdef WSSMP_DBG
1458        chk(32, j2, adjln);
1459        chk(33, j0, adjln);
1460#endif
1461        adjncy[i4 - 1] = adjncy[j2 - 1];
1462        adjncy[j2 - 1] = adjncy[j0 - 1];
1463        adjncy[j0 - 1] = node;
1464        lsize[i - 1] = i4 - j0 + 1;
1465        i5++;
1466#ifdef WSSMP_DBG
1467        chk(35, i5, n);
1468#endif
1469        j = head[i5 - 1];
1470        if (j > 0) {
1471          snxt[i - 1] = perm[j - 1];
1472          perm[j - 1] = i;
1473        } else {
1474          snxt[i - 1] = -j;
1475          head[i5 - 1] = -i;
1476        }
1477        perm[i - 1] = i5;
1478      }
1479    }
1480#ifdef WSSMP_DBG
1481    chk(36, node, n);
1482#endif
1483    dgree[node - 1] = nodeg;
1484    if (maxmum < nodeg)
1485      maxmum = nodeg;
1486    fltag += maxmum;
1487#ifdef WSSMP_DBG
1488    if (fltag + n + 1 < 0)
1489      printf("ERR**: fltag = %d (A)\n", fltag);
1490#endif
1491    for (j3 = j4; j3 < j5; ++j3) {
1492#ifdef WSSMP_DBG
1493      chk(37, j3, adjln);
1494#endif
1495      i = adjncy[j3 - 1];
1496#ifdef WSSMP_DBG
1497      chk(38, i, n);
1498#endif
1499      if (varbl[i - 1] < 0) {
1500        i5 = perm[i - 1];
1501#ifdef WSSMP_DBG
1502        chk(39, i5, n);
1503#endif
1504        j = head[i5 - 1];
1505        if (j) {
1506          if (j < 0) {
1507            head[i5 - 1] = 0, i = -j;
1508          } else {
1509#ifdef WSSMP_DBG
1510            chk(40, j, n);
1511#endif
1512            i = perm[j - 1];
1513            perm[j - 1] = 0;
1514          }
1515          while (i) {
1516#ifdef WSSMP_DBG
1517            chk(41, i, n);
1518#endif
1519            if (!snxt[i - 1]) {
1520              i = 0;
1521            } else {
1522              k = invp[i - 1];
1523              i2 = xadj[i - 1] + (scln = lsize[i - 1]);
1524              for (l = xadj[i - 1] + 1; l < i2; l++) {
1525#ifdef WSSMP_DBG
1526                chk(42, l, adjln);
1527                chk(43, adjncy[l - 1], n);
1528#endif
1529                flag[adjncy[l - 1] - 1] = fltag;
1530              }
1531              jpp = i;
1532              j = snxt[i - 1];
1533              while (j) {
1534#ifdef WSSMP_DBG
1535                chk(44, j, n);
1536#endif
1537                if (lsize[j - 1] == scln && invp[j - 1] == k) {
1538                  i2 = xadj[j - 1] + scln;
1539                  j8 = 1;
1540                  for (l = xadj[j - 1] + 1; l < i2; l++) {
1541#ifdef WSSMP_DBG
1542                    chk(45, l, adjln);
1543                    chk(46, adjncy[l - 1], n);
1544#endif
1545                    if (flag[adjncy[l - 1] - 1] != fltag) {
1546                      j8 = 0;
1547                      break;
1548                    }
1549                  }
1550                  if (j8) {
1551                    xadj[j - 1] = -i;
1552                    varbl[i - 1] += varbl[j - 1];
1553                    varbl[j - 1] = invp[j - 1] = 0;
1554#ifdef WSSMP_DBG
1555                    chk(47, j, n);
1556                    chk(48, jpp, n);
1557#endif
1558                    j = snxt[j - 1];
1559                    snxt[jpp - 1] = j;
1560                  } else {
1561                    jpp = j;
1562#ifdef WSSMP_DBG
1563                    chk(49, j, n);
1564#endif
1565                    j = snxt[j - 1];
1566                  }
1567                } else {
1568                  jpp = j;
1569#ifdef WSSMP_DBG
1570                  chk(50, j, n);
1571#endif
1572                  j = snxt[j - 1];
1573                }
1574              }
1575              fltag++;
1576#ifdef WSSMP_DBG
1577              if (fltag + n + 1 < 0)
1578                printf("ERR**: fltag = %d (B)\n", fltag);
1579#endif
1580#ifdef WSSMP_DBG
1581              chk(51, i, n);
1582#endif
1583              i = snxt[i - 1];
1584            }
1585          }
1586        }
1587      }
1588    }
1589    j8 = nm1 - counter;
1590    switch (speed) {
1591    case 1:
1592      for (j = j3 = j4; j3 < j5; j3++) {
1593#ifdef WSSMP_DBG
1594        chk(52, j3, adjln);
1595#endif
1596        i = adjncy[j3 - 1];
1597#ifdef WSSMP_DBG
1598        chk(53, i, n);
1599#endif
1600        numii = varbl[i - 1];
1601        if (numii < 0) {
1602          k = 0;
1603          i4 = (l = xadj[i - 1]) + invp[i - 1];
1604          for (; l < i4; l++) {
1605#ifdef WSSMP_DBG
1606            chk(54, l, adjln);
1607            chk(55, adjncy[l - 1], n);
1608#endif
1609            i5 = dgree[adjncy[l - 1] - 1];
1610            if (k < i5)
1611              k = i5;
1612          }
1613          x = static_cast< double >(k - 1);
1614          varbl[i - 1] = abs(numii);
1615          j9 = dgree[i - 1] + nodeg;
1616          deg = (j8 > j9 ? j9 : j8) + numii;
1617          if (deg < 1)
1618            deg = 1;
1619          y = static_cast< double >(deg);
1620          dgree[i - 1] = deg;
1621          /*
1622                                    printf("%i %f should >= %i %f\n",deg,y,k-1,x);
1623                                    if (y < x) printf("** \n"); else printf("\n");
1624                         */
1625          y = y * y - y;
1626          x = y - (x * x - x);
1627          if (x < 1.1)
1628            x = 1.1;
1629          deg = static_cast< WSI >(sqrt(x));
1630          /*         if (deg < 1) deg = 1; */
1631          erscore[i - 1] = deg;
1632#ifdef WSSMP_DBG
1633          chk(56, deg, n - 1);
1634#endif
1635          nnode = head[deg - 1];
1636          if (nnode)
1637            perm[nnode - 1] = i;
1638          snxt[i - 1] = nnode;
1639          perm[i - 1] = 0;
1640#ifdef WSSMP_DBG
1641          chk(57, j, adjln);
1642#endif
1643          head[deg - 1] = adjncy[j - 1] = i;
1644          j++;
1645          if (deg < mindeg)
1646            mindeg = deg;
1647        }
1648      }
1649      break;
1650    case 2:
1651      for (j = j3 = j4; j3 < j5; j3++) {
1652#ifdef WSSMP_DBG
1653        chk(58, j3, adjln);
1654#endif
1655        i = adjncy[j3 - 1];
1656#ifdef WSSMP_DBG
1657        chk(59, i, n);
1658#endif
1659        numii = varbl[i - 1];
1660        if (numii < 0) {
1661          x = static_cast< double >(dgree[adjncy[xadj[i - 1] - 1] - 1] - 1);
1662          varbl[i - 1] = abs(numii);
1663          j9 = dgree[i - 1] + nodeg;
1664          deg = (j8 > j9 ? j9 : j8) + numii;
1665          if (deg < 1)
1666            deg = 1;
1667          y = static_cast< double >(deg);
1668          dgree[i - 1] = deg;
1669          /*
1670                                    printf("%i %f should >= %i %f",deg,y,dgree[adjncy[xadj[i - 1] - 1] - 1]-1,x);
1671                                    if (y < x) printf("** \n"); else printf("\n");
1672                         */
1673          y = y * y - y;
1674          x = y - (x * x - x);
1675          if (x < 1.1)
1676            x = 1.1;
1677          deg = static_cast< WSI >(sqrt(x));
1678          /*         if (deg < 1) deg = 1; */
1679          erscore[i - 1] = deg;
1680#ifdef WSSMP_DBG
1681          chk(60, deg, n - 1);
1682#endif
1683          nnode = head[deg - 1];
1684          if (nnode)
1685            perm[nnode - 1] = i;
1686          snxt[i - 1] = nnode;
1687          perm[i - 1] = 0;
1688#ifdef WSSMP_DBG
1689          chk(61, j, adjln);
1690#endif
1691          head[deg - 1] = adjncy[j - 1] = i;
1692          j++;
1693          if (deg < mindeg)
1694            mindeg = deg;
1695        }
1696      }
1697      break;
1698    default:
1699      for (j = j3 = j4; j3 < j5; j3++) {
1700#ifdef WSSMP_DBG
1701        chk(62, j3, adjln);
1702#endif
1703        i = adjncy[j3 - 1];
1704#ifdef WSSMP_DBG
1705        chk(63, i, n);
1706#endif
1707        numii = varbl[i - 1];
1708        if (numii < 0) {
1709          varbl[i - 1] = abs(numii);
1710          j9 = dgree[i - 1] + nodeg;
1711          deg = (j8 > j9 ? j9 : j8) + numii;
1712          if (deg < 1)
1713            deg = 1;
1714          dgree[i - 1] = deg;
1715#ifdef WSSMP_DBG
1716          chk(64, deg, n - 1);
1717#endif
1718          nnode = head[deg - 1];
1719          if (nnode)
1720            perm[nnode - 1] = i;
1721          snxt[i - 1] = nnode;
1722          perm[i - 1] = 0;
1723#ifdef WSSMP_DBG
1724          chk(65, j, adjln);
1725#endif
1726          head[deg - 1] = adjncy[j - 1] = i;
1727          j++;
1728          if (deg < mindeg)
1729            mindeg = deg;
1730        }
1731      }
1732      break;
1733    }
1734    if (currloc) {
1735#ifdef WSSMP_DBG
1736      chk(66, node, n);
1737#endif
1738      locatns += (lsize[node - 1] - currloc), locaux = j;
1739    }
1740    varbl[node - 1] = numpiv + nodeg;
1741    lsize[node - 1] = j - j4;
1742    if (!lsize[node - 1])
1743      flag[node - 1] = xadj[node - 1] = 0;
1744  }
1745  for (l = 1; l <= n; l++) {
1746    if (!invp[l - 1]) {
1747      for (i = -xadj[l - 1]; invp[i - 1] >= 0; i = -xadj[i - 1]) {
1748      };
1749      tn = i;
1750#ifdef WSSMP_DBG
1751      chk(67, tn, n);
1752#endif
1753      k = -invp[tn - 1];
1754      i = l;
1755#ifdef WSSMP_DBG
1756      chk(68, i, n);
1757#endif
1758      while (invp[i - 1] >= 0) {
1759        nnode = -xadj[i - 1];
1760        xadj[i - 1] = -tn;
1761        if (!invp[i - 1])
1762          invp[i - 1] = k++;
1763        i = nnode;
1764      }
1765      invp[tn - 1] = -k;
1766    }
1767  }
1768  for (l = 0; l < n; l++) {
1769    i = abs(invp[l]);
1770#ifdef WSSMP_DBG
1771    chk(69, i, n);
1772#endif
1773    invp[l] = i;
1774    perm[i - 1] = l + 1;
1775  }
1776  return;
1777}
1778// This code is not needed, but left in in case needed sometime
1779#if 0
1780/*C--------------------------------------------------------------------------*/
1781void amlfdr(WSI *n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI *adjln,
1782            WSI *locaux, WSI varbl[], WSI snxt[], WSI perm[],
1783            WSI head[], WSI invp[], WSI lsize[], WSI flag[], WSI *ispeed)
1784{
1785     WSI nn, nlocaux, nadjln, speed, i, j, mx, mxj, *erscore;
1786
1787#ifdef WSSMP_DBG
1788     printf("Calling amlfdr for n, speed = %d, %d\n", *n, *ispeed);
1789#endif
1790
1791     if ((nn = *n) == 0) return;
1792
1793#ifdef WSSMP_DBG
1794     if (nn == 31) {
1795          printf("n = %d; adjln = %d; locaux = %d; ispeed = %d\n",
1796                 *n, *adjln, *locaux, *ispeed);
1797     }
1798#endif
1799
1800     if (nn < NORDTHRESH) {
1801          for (i = 0; i < nn; i++) lsize[i] = i;
1802          for (i = nn; i > 0; i--) {
1803               mx = dgree[0];
1804               mxj = 0;
1805               for (j = 1; j < i; j++)
1806                    if (dgree[j] > mx) {
1807                         mx = dgree[j];
1808                         mxj = j;
1809                    }
1810               invp[lsize[mxj]] = i;
1811               dgree[mxj] = dgree[i-1];
1812               lsize[mxj] = lsize[i-1];
1813          }
1814          return;
1815     }
1816     speed = *ispeed;
1817     if (speed < 3) {
1818          /*
1819              erscore = (WSI *)malloc(nn * sizeof(WSI));
1820              if (erscore == NULL) speed = 3;
1821          */
1822          wscmal (&nn, &i, &erscore);
1823          if (i != 0) speed = 3;
1824     }
1825     if (speed > 2) erscore = dgree;
1826     if (speed < 3) {
1827          for (i = 0; i < nn; i++) {
1828               perm[i] = 0;
1829               invp[i] = 0;
1830               head[i] = 0;
1831               flag[i] = 1;
1832               varbl[i] = 1;
1833               lsize[i] = dgree[i];
1834               erscore[i] = dgree[i];
1835          }
1836     } else {
1837          for (i = 0; i < nn; i++) {
1838               perm[i] = 0;
1839               invp[i] = 0;
1840               head[i] = 0;
1841               flag[i] = 1;
1842               varbl[i] = 1;
1843               lsize[i] = dgree[i];
1844          }
1845     }
1846     nlocaux = *locaux;
1847     nadjln = *adjln;
1848
1849     myamlf(nn, xadj, adjncy, dgree, varbl, snxt, perm, invp,
1850            head, lsize, flag, erscore, nlocaux, nadjln, speed);
1851     /*
1852       if (speed < 3) free(erscore);
1853     */
1854     if (speed < 3) wscfr(&erscore);
1855     return;
1856}
1857#endif // end of taking out amlfdr
1858/*C--------------------------------------------------------------------------*/
1859#endif
1860// Orders rows
1861int ClpCholeskyBase::orderAMD()
1862{
1863  permuteInverse_ = new int[numberRows_];
1864  permute_ = new int[numberRows_];
1865  // Do ordering
1866  int returnCode = 0;
1867  // get full matrix
1868  int space = 2 * sizeFactor_ + 10000 + 4 * numberRows_;
1869  int *temp = new int[space];
1870  int *count = new int[numberRows_];
1871  int *tempStart = new int[numberRows_ + 1];
1872  memset(count, 0, numberRows_ * sizeof(int));
1873  for (int iRow = 0; iRow < numberRows_; iRow++) {
1874    count[iRow] += static_cast< int >(choleskyStart_[iRow + 1] - choleskyStart_[iRow] - 1);
1875    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
1876      int jRow = choleskyRow_[j];
1877      count[jRow]++;
1878    }
1879  }
1880#define OFFSET 1
1881  int sizeFactor = 0;
1882  for (int iRow = 0; iRow < numberRows_; iRow++) {
1883    int length = count[iRow];
1884    permute_[iRow] = length;
1885    // add 1 to starts
1886    tempStart[iRow] = sizeFactor + OFFSET;
1887    count[iRow] = sizeFactor;
1888    sizeFactor += length;
1889  }
1890  tempStart[numberRows_] = sizeFactor + OFFSET;
1891  // add 1 to rows
1892  for (int iRow = 0; iRow < numberRows_; iRow++) {
1893    assert(choleskyRow_[choleskyStart_[iRow]] == iRow);
1894    for (CoinBigIndex j = choleskyStart_[iRow] + 1; j < choleskyStart_[iRow + 1]; j++) {
1895      int jRow = choleskyRow_[j];
1896      int put = count[iRow];
1897      temp[put] = jRow + OFFSET;
1898      count[iRow]++;
1899      put = count[jRow];
1900      temp[put] = iRow + OFFSET;
1901      count[jRow]++;
1902    }
1903  }
1904  for (int iRow = 1; iRow < numberRows_; iRow++)
1905    assert(count[iRow - 1] == tempStart[iRow] - OFFSET);
1906  delete[] choleskyRow_;
1907  choleskyRow_ = temp;
1908  delete[] choleskyStart_;
1909  choleskyStart_ = tempStart;
1910  int locaux = sizeFactor + OFFSET;
1911  delete[] count;
1912  int speed = integerParameters_[0];
1913  if (speed < 1 || speed > 2)
1914    speed = 3;
1915  int *use = new int[((speed < 3) ? 7 : 6) * numberRows_];
1916  int *dgree = use;
1917  int *varbl = dgree + numberRows_;
1918  int *snxt = varbl + numberRows_;
1919  int *head = snxt + numberRows_;
1920  int *lsize = head + numberRows_;
1921  int *flag = lsize + numberRows_;
1922  int *erscore;
1923  for (int i = 0; i < numberRows_; i++) {
1924    dgree[i] = choleskyStart_[i + 1] - choleskyStart_[i];
1925    head[i] = dgree[i];
1926    snxt[i] = 0;
1927    permute_[i] = 0;
1928    permuteInverse_[i] = 0;
1929    head[i] = 0;
1930    flag[i] = 1;
1931    varbl[i] = 1;
1932    lsize[i] = dgree[i];
1933  }
1934  if (speed < 3) {
1935    erscore = flag + numberRows_;
1936    for (int i = 0; i < numberRows_; i++)
1937      erscore[i] = dgree[i];
1938  } else {
1939    erscore = dgree;
1940  }
1941  myamlf(numberRows_, choleskyStart_, choleskyRow_,
1942    dgree, varbl, snxt, permute_, permuteInverse_,
1943    head, lsize, flag, erscore, locaux, space, speed);
1944  for (int iRow = 0; iRow < numberRows_; iRow++) {
1945    permute_[iRow]--;
1946  }
1947  for (int iRow = 0; iRow < numberRows_; iRow++) {
1948    permuteInverse_[permute_[iRow]] = iRow;
1949  }
1950  for (int iRow = 0; iRow < numberRows_; iRow++) {
1951    assert(permuteInverse_[iRow] >= 0 && permuteInverse_[iRow] < numberRows_);
1952  }
1953  delete[] use;
1954  delete[] choleskyRow_;
1955  choleskyRow_ = NULL;
1956  delete[] choleskyStart_;
1957  choleskyStart_ = NULL;
1958  return returnCode;
1959}
1960/* Does Symbolic factorization given permutation.
1961   This is called immediately after order.  If user provides this then
1962   user must provide factorize and solve.  Otherwise the default factorization is used
1963   returns non-zero if not enough memory */
1964int ClpCholeskyBase::symbolic()
1965{
1966  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
1967  const int *columnLength = model_->clpMatrix()->getVectorLengths();
1968  const int *row = model_->clpMatrix()->getIndices();
1969  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
1970  const int *rowLength = rowCopy_->getVectorLengths();
1971  const int *column = rowCopy_->getIndices();
1972  int numberRowsModel = model_->numberRows();
1973  int numberColumns = model_->numberColumns();
1974  int numberTotal = numberColumns + numberRowsModel;
1975  CoinPackedMatrix *quadratic = NULL;
1976  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
1977  if (quadraticObj)
1978    quadratic = quadraticObj->quadraticObjective();
1979  // We need an array for counts
1980  int *used = new int[numberRows_ + 1];
1981  // If KKT then re-order so negative first
1982  if (doKKT_) {
1983    int nn = 0;
1984    int np = 0;
1985    int iRow;
1986    for (iRow = 0; iRow < numberRows_; iRow++) {
1987      int originalRow = permute_[iRow];
1988      if (originalRow < numberTotal)
1989        permute_[nn++] = originalRow;
1990      else
1991        used[np++] = originalRow;
1992    }
1993    CoinMemcpyN(used, np, permute_ + nn);
1994    for (iRow = 0; iRow < numberRows_; iRow++)
1995      permuteInverse_[permute_[iRow]] = iRow;
1996  }
1997  CoinZeroN(used, numberRows_);
1998  int iRow;
1999  int iColumn;
2000  bool noMemory = false;
2001  int *Astart = new int[numberRows_ + 1];
2002  int *Arow = NULL;
2003  try {
2004    Arow = new int[sizeFactor_];
2005  } catch (...) {
2006    // no memory
2007    delete[] Astart;
2008    return -1;
2009  }
2010  choleskyStart_ = new int[numberRows_ + 1];
2011  link_ = new int[numberRows_];
2012  workInteger_ = new int[numberRows_];
2013  indexStart_ = new int[numberRows_];
2014  clique_ = new int[numberRows_];
2015  // Redo so permuted upper triangular
2016  sizeFactor_ = 0;
2017  int *which = Arow;
2018  if (!doKKT_) {
2019    for (iRow = 0; iRow < numberRows_; iRow++) {
2020      int number = 0;
2021      int iOriginalRow = permute_[iRow];
2022      Astart[iRow] = sizeFactor_;
2023      CoinBigIndex startRow = rowStart[iOriginalRow];
2024      CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
2025      for (CoinBigIndex k = startRow; k < endRow; k++) {
2026        int iColumn = column[k];
2027        if (!whichDense_ || !whichDense_[iColumn]) {
2028          CoinBigIndex start = columnStart[iColumn];
2029          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2030          for (CoinBigIndex j = start; j < end; j++) {
2031            int jRow = row[j];
2032            int jNewRow = permuteInverse_[jRow];
2033            if (jNewRow < iRow) {
2034              if (!used[jNewRow]) {
2035                used[jNewRow] = 1;
2036                which[number++] = jNewRow;
2037              }
2038            }
2039          }
2040        }
2041      }
2042      sizeFactor_ += number;
2043      int j;
2044      for (j = 0; j < number; j++)
2045        used[which[j]] = 0;
2046      // Sort
2047      std::sort(which, which + number);
2048      // move which on
2049      which += number;
2050    }
2051  } else {
2052    // KKT
2053    // transpose
2054    ClpMatrixBase *rowCopy = model_->clpMatrix()->reverseOrderedCopy();
2055    const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
2056    const int *rowLength = rowCopy->getVectorLengths();
2057    const int *column = rowCopy->getIndices();
2058    // temp
2059    bool permuted = false;
2060    for (iRow = 0; iRow < numberRows_; iRow++) {
2061      if (permute_[iRow] != iRow) {
2062        permuted = true;
2063        break;
2064      }
2065    }
2066    if (permuted) {
2067      // Need to permute - ugly
2068      if (!quadratic) {
2069        for (iRow = 0; iRow < numberRows_; iRow++) {
2070          Astart[iRow] = sizeFactor_;
2071          int iOriginalRow = permute_[iRow];
2072          if (iOriginalRow < numberColumns) {
2073            // A may be upper triangular by mistake
2074            iColumn = iOriginalRow;
2075            CoinBigIndex start = columnStart[iColumn];
2076            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2077            for (CoinBigIndex j = start; j < end; j++) {
2078              int kRow = row[j] + numberTotal;
2079              kRow = permuteInverse_[kRow];
2080              if (kRow < iRow)
2081                Arow[sizeFactor_++] = kRow;
2082            }
2083          } else if (iOriginalRow < numberTotal) {
2084            int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
2085            if (kRow < iRow)
2086              Arow[sizeFactor_++] = kRow;
2087          } else {
2088            int kRow = iOriginalRow - numberTotal;
2089            CoinBigIndex start = rowStart[kRow];
2090            CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
2091            for (CoinBigIndex j = start; j < end; j++) {
2092              int jRow = column[j];
2093              int jNewRow = permuteInverse_[jRow];
2094              if (jNewRow < iRow)
2095                Arow[sizeFactor_++] = jNewRow;
2096            }
2097            // slack - should it be permute
2098            kRow = permuteInverse_[kRow + numberColumns];
2099            if (kRow < iRow)
2100              Arow[sizeFactor_++] = kRow;
2101          }
2102          // Sort
2103          std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
2104        }
2105      } else {
2106        // quadratic
2107        // transpose
2108        CoinPackedMatrix quadraticT;
2109        quadraticT.reverseOrderedCopyOf(*quadratic);
2110        const int *columnQuadratic = quadraticT.getIndices();
2111        const CoinBigIndex *columnQuadraticStart = quadraticT.getVectorStarts();
2112        const int *columnQuadraticLength = quadraticT.getVectorLengths();
2113        for (iRow = 0; iRow < numberRows_; iRow++) {
2114          Astart[iRow] = sizeFactor_;
2115          int iOriginalRow = permute_[iRow];
2116          if (iOriginalRow < numberColumns) {
2117            // Quadratic bit
2118            CoinBigIndex j;
2119            for (j = columnQuadraticStart[iOriginalRow];
2120                 j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
2121              int jColumn = columnQuadratic[j];
2122              int jNewColumn = permuteInverse_[jColumn];
2123              if (jNewColumn < iRow)
2124                Arow[sizeFactor_++] = jNewColumn;
2125            }
2126            // A may be upper triangular by mistake
2127            iColumn = iOriginalRow;
2128            CoinBigIndex start = columnStart[iColumn];
2129            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2130            for (j = start; j < end; j++) {
2131              int kRow = row[j] + numberTotal;
2132              kRow = permuteInverse_[kRow];
2133              if (kRow < iRow)
2134                Arow[sizeFactor_++] = kRow;
2135            }
2136          } else if (iOriginalRow < numberTotal) {
2137            int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
2138            if (kRow < iRow)
2139              Arow[sizeFactor_++] = kRow;
2140          } else {
2141            int kRow = iOriginalRow - numberTotal;
2142            CoinBigIndex start = rowStart[kRow];
2143            CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
2144            for (CoinBigIndex j = start; j < end; j++) {
2145              int jRow = column[j];
2146              int jNewRow = permuteInverse_[jRow];
2147              if (jNewRow < iRow)
2148                Arow[sizeFactor_++] = jNewRow;
2149            }
2150            // slack - should it be permute
2151            kRow = permuteInverse_[kRow + numberColumns];
2152            if (kRow < iRow)
2153              Arow[sizeFactor_++] = kRow;
2154          }
2155          // Sort
2156          std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
2157        }
2158      }
2159    } else {
2160      if (!quadratic) {
2161        for (iRow = 0; iRow < numberRows_; iRow++) {
2162          Astart[iRow] = sizeFactor_;
2163        }
2164      } else {
2165        // Quadratic
2166        // transpose
2167        CoinPackedMatrix quadraticT;
2168        quadraticT.reverseOrderedCopyOf(*quadratic);
2169        const int *columnQuadratic = quadraticT.getIndices();
2170        const CoinBigIndex *columnQuadraticStart = quadraticT.getVectorStarts();
2171        const int *columnQuadraticLength = quadraticT.getVectorLengths();
2172        //const double * quadraticElement = quadraticT.getElements();
2173        for (iRow = 0; iRow < numberColumns; iRow++) {
2174          int iOriginalRow = permute_[iRow];
2175          Astart[iRow] = sizeFactor_;
2176          for (CoinBigIndex j = columnQuadraticStart[iOriginalRow];
2177               j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
2178            int jColumn = columnQuadratic[j];
2179            int jNewColumn = permuteInverse_[jColumn];
2180            if (jNewColumn < iRow)
2181              Arow[sizeFactor_++] = jNewColumn;
2182          }
2183        }
2184      }
2185      int iRow;
2186      // slacks
2187      for (iRow = 0; iRow < numberRowsModel; iRow++) {
2188        Astart[iRow + numberColumns] = sizeFactor_;
2189      }
2190      for (iRow = 0; iRow < numberRowsModel; iRow++) {
2191        Astart[iRow + numberTotal] = sizeFactor_;
2192        CoinBigIndex start = rowStart[iRow];
2193        CoinBigIndex end = rowStart[iRow] + rowLength[iRow];
2194        for (CoinBigIndex j = start; j < end; j++) {
2195          Arow[sizeFactor_++] = column[j];
2196        }
2197        // slack
2198        Arow[sizeFactor_++] = numberColumns + iRow;
2199      }
2200    }
2201    delete rowCopy;
2202  }
2203  Astart[numberRows_] = sizeFactor_;
2204  firstDense_ = numberRows_;
2205  symbolic1(Astart, Arow);
2206  // Now fill in indices
2207  try {
2208    // too big
2209    choleskyRow_ = new int[sizeFactor_];
2210  } catch (...) {
2211    // no memory
2212    noMemory = true;
2213  }
2214  double sizeFactor = sizeFactor_;
2215  if (!noMemory) {
2216    // Do lower triangular
2217    sizeFactor_ = 0;
2218    int *which = Arow;
2219    if (!doKKT_) {
2220      for (iRow = 0; iRow < numberRows_; iRow++) {
2221        int number = 0;
2222        int iOriginalRow = permute_[iRow];
2223        Astart[iRow] = sizeFactor_;
2224        if (!rowsDropped_[iOriginalRow]) {
2225          CoinBigIndex startRow = rowStart[iOriginalRow];
2226          CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
2227          for (CoinBigIndex k = startRow; k < endRow; k++) {
2228            int iColumn = column[k];
2229            if (!whichDense_ || !whichDense_[iColumn]) {
2230              CoinBigIndex start = columnStart[iColumn];
2231              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2232              for (CoinBigIndex j = start; j < end; j++) {
2233                int jRow = row[j];
2234                int jNewRow = permuteInverse_[jRow];
2235                if (jNewRow > iRow && !rowsDropped_[jRow]) {
2236                  if (!used[jNewRow]) {
2237                    used[jNewRow] = 1;
2238                    which[number++] = jNewRow;
2239                  }
2240                }
2241              }
2242            }
2243          }
2244          sizeFactor_ += number;
2245          int j;
2246          for (j = 0; j < number; j++)
2247            used[which[j]] = 0;
2248          // Sort
2249          std::sort(which, which + number);
2250          // move which on
2251          which += number;
2252        }
2253      }
2254    } else {
2255      // KKT
2256      // temp
2257      bool permuted = false;
2258      for (iRow = 0; iRow < numberRows_; iRow++) {
2259        if (permute_[iRow] != iRow) {
2260          permuted = true;
2261          break;
2262        }
2263      }
2264      // but fake it
2265      for (iRow = 0; iRow < numberRows_; iRow++) {
2266        //permute_[iRow]=iRow; // force no permute
2267        //permuteInverse_[permute_[iRow]]=iRow;
2268      }
2269      if (permuted) {
2270        // Need to permute - ugly
2271        if (!quadratic) {
2272          for (iRow = 0; iRow < numberRows_; iRow++) {
2273            Astart[iRow] = sizeFactor_;
2274            int iOriginalRow = permute_[iRow];
2275            if (iOriginalRow < numberColumns) {
2276              iColumn = iOriginalRow;
2277              CoinBigIndex start = columnStart[iColumn];
2278              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2279              for (CoinBigIndex j = start; j < end; j++) {
2280                int kRow = row[j] + numberTotal;
2281                kRow = permuteInverse_[kRow];
2282                if (kRow > iRow)
2283                  Arow[sizeFactor_++] = kRow;
2284              }
2285            } else if (iOriginalRow < numberTotal) {
2286              int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
2287              if (kRow > iRow)
2288                Arow[sizeFactor_++] = kRow;
2289            } else {
2290              int kRow = iOriginalRow - numberTotal;
2291              CoinBigIndex start = rowStart[kRow];
2292              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
2293              for (CoinBigIndex j = start; j < end; j++) {
2294                int jRow = column[j];
2295                int jNewRow = permuteInverse_[jRow];
2296                if (jNewRow > iRow)
2297                  Arow[sizeFactor_++] = jNewRow;
2298              }
2299              // slack - should it be permute
2300              kRow = permuteInverse_[kRow + numberColumns];
2301              if (kRow > iRow)
2302                Arow[sizeFactor_++] = kRow;
2303            }
2304            // Sort
2305            std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
2306          }
2307        } else {
2308          // quadratic
2309          const int *columnQuadratic = quadratic->getIndices();
2310          const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
2311          const int *columnQuadraticLength = quadratic->getVectorLengths();
2312          for (iRow = 0; iRow < numberRows_; iRow++) {
2313            Astart[iRow] = sizeFactor_;
2314            int iOriginalRow = permute_[iRow];
2315            if (iOriginalRow < numberColumns) {
2316              // Quadratic bit
2317              CoinBigIndex j;
2318              for (j = columnQuadraticStart[iOriginalRow];
2319                   j < columnQuadraticStart[iOriginalRow] + columnQuadraticLength[iOriginalRow]; j++) {
2320                int jColumn = columnQuadratic[j];
2321                int jNewColumn = permuteInverse_[jColumn];
2322                if (jNewColumn > iRow)
2323                  Arow[sizeFactor_++] = jNewColumn;
2324              }
2325              iColumn = iOriginalRow;
2326              CoinBigIndex start = columnStart[iColumn];
2327              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2328              for (j = start; j < end; j++) {
2329                int kRow = row[j] + numberTotal;
2330                kRow = permuteInverse_[kRow];
2331                if (kRow > iRow)
2332                  Arow[sizeFactor_++] = kRow;
2333              }
2334            } else if (iOriginalRow < numberTotal) {
2335              int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
2336              if (kRow > iRow)
2337                Arow[sizeFactor_++] = kRow;
2338            } else {
2339              int kRow = iOriginalRow - numberTotal;
2340              CoinBigIndex start = rowStart[kRow];
2341              CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
2342              for (CoinBigIndex j = start; j < end; j++) {
2343                int jRow = column[j];
2344                int jNewRow = permuteInverse_[jRow];
2345                if (jNewRow > iRow)
2346                  Arow[sizeFactor_++] = jNewRow;
2347              }
2348              // slack - should it be permute
2349              kRow = permuteInverse_[kRow + numberColumns];
2350              if (kRow > iRow)
2351                Arow[sizeFactor_++] = kRow;
2352            }
2353            // Sort
2354            std::sort(Arow + Astart[iRow], Arow + sizeFactor_);
2355          }
2356        }
2357      } else {
2358        if (!quadratic) {
2359          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2360            Astart[iColumn] = sizeFactor_;
2361            CoinBigIndex start = columnStart[iColumn];
2362            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2363            for (CoinBigIndex j = start; j < end; j++) {
2364              Arow[sizeFactor_++] = row[j] + numberTotal;
2365            }
2366          }
2367        } else {
2368          // Quadratic
2369          const int *columnQuadratic = quadratic->getIndices();
2370          const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
2371          const int *columnQuadraticLength = quadratic->getVectorLengths();
2372          //const double * quadraticElement = quadratic->getElements();
2373          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2374            Astart[iColumn] = sizeFactor_;
2375            CoinBigIndex j;
2376            for (j = columnQuadraticStart[iColumn];
2377                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
2378              int jColumn = columnQuadratic[j];
2379              if (jColumn > iColumn)
2380                Arow[sizeFactor_++] = jColumn;
2381            }
2382            CoinBigIndex start = columnStart[iColumn];
2383            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2384            for (j = start; j < end; j++) {
2385              Arow[sizeFactor_++] = row[j] + numberTotal;
2386            }
2387          }
2388        }
2389        // slacks
2390        for (iRow = 0; iRow < numberRowsModel; iRow++) {
2391          Astart[iRow + numberColumns] = sizeFactor_;
2392          Arow[sizeFactor_++] = iRow + numberTotal;
2393        }
2394        // Transpose - nonzero diagonal (may regularize)
2395        for (iRow = 0; iRow < numberRowsModel; iRow++) {
2396          Astart[iRow + numberTotal] = sizeFactor_;
2397        }
2398      }
2399      Astart[numberRows_] = sizeFactor_;
2400    }
2401    symbolic2(Astart, Arow);
2402    if (sizeIndex_ < sizeFactor_) {
2403      int *indices = NULL;
2404      try {
2405        indices = new int[sizeIndex_];
2406      } catch (...) {
2407        // no memory
2408        noMemory = true;
2409      }
2410      if (!noMemory) {
2411        CoinMemcpyN(choleskyRow_, sizeIndex_, indices);
2412        delete[] choleskyRow_;
2413        choleskyRow_ = indices;
2414      }
2415    }
2416  }
2417  delete[] used;
2418  // Use cholesky regions
2419  delete[] Astart;
2420  delete[] Arow;
2421  double flops = 0.0;
2422  for (iRow = 0; iRow < numberRows_; iRow++) {
2423    int length = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
2424    flops += static_cast< double >(length) * (length + 2.0);
2425  }
2426  if (model_->messageHandler()->logLevel() > 0)
2427    std::cout << sizeFactor << " elements in sparse Cholesky, flop count " << flops << std::endl;
2428  try {
2429    sparseFactor_ = new longDouble[sizeFactor_];
2430#if CLP_LONG_CHOLESKY != 1
2431    workDouble_ = new longDouble[numberRows_];
2432#else
2433    // actually long double
2434    workDouble_ = reinterpret_cast< double * >(new CoinWorkDouble[numberRows_]);
2435#endif
2436    diagonal_ = new longDouble[numberRows_];
2437  } catch (...) {
2438    // no memory
2439    noMemory = true;
2440  }
2441  if (noMemory) {
2442    delete[] choleskyRow_;
2443    choleskyRow_ = NULL;
2444    delete[] choleskyStart_;
2445    choleskyStart_ = NULL;
2446    delete[] permuteInverse_;
2447    permuteInverse_ = NULL;
2448    delete[] permute_;
2449    permute_ = NULL;
2450    delete[] choleskyStart_;
2451    choleskyStart_ = NULL;
2452    delete[] indexStart_;
2453    indexStart_ = NULL;
2454    delete[] link_;
2455    link_ = NULL;
2456    delete[] workInteger_;
2457    workInteger_ = NULL;
2458    delete[] sparseFactor_;
2459    sparseFactor_ = NULL;
2460    delete[] workDouble_;
2461    workDouble_ = NULL;
2462    delete[] diagonal_;
2463    diagonal_ = NULL;
2464    delete[] clique_;
2465    clique_ = NULL;
2466    return -1;
2467  }
2468  return 0;
2469}
2470int ClpCholeskyBase::symbolic1(const int *Astart, const int *Arow)
2471{
2472  int *marked = reinterpret_cast< int * >(workInteger_);
2473  int iRow;
2474  // may not need to do this here but makes debugging easier
2475  for (iRow = 0; iRow < numberRows_; iRow++) {
2476    marked[iRow] = -1;
2477    link_[iRow] = -1;
2478    choleskyStart_[iRow] = 0; // counts
2479  }
2480  for (iRow = 0; iRow < numberRows_; iRow++) {
2481    marked[iRow] = iRow;
2482    for (CoinBigIndex j = Astart[iRow]; j < Astart[iRow + 1]; j++) {
2483      int kRow = Arow[j];
2484      while (marked[kRow] != iRow) {
2485        if (link_[kRow] < 0)
2486          link_[kRow] = iRow;
2487        choleskyStart_[kRow]++;
2488        marked[kRow] = iRow;
2489        kRow = link_[kRow];
2490      }
2491    }
2492  }
2493  sizeFactor_ = 0;
2494  for (iRow = 0; iRow < numberRows_; iRow++) {
2495    int number = choleskyStart_[iRow];
2496    choleskyStart_[iRow] = sizeFactor_;
2497    sizeFactor_ += number;
2498  }
2499  choleskyStart_[numberRows_] = sizeFactor_;
2500  return sizeFactor_;
2501  ;
2502}
2503void ClpCholeskyBase::symbolic2(const int *Astart, const int *Arow)
2504{
2505  int *mergeLink = clique_;
2506  int *marker = reinterpret_cast< int * >(workInteger_);
2507  int iRow;
2508  for (iRow = 0; iRow < numberRows_; iRow++) {
2509    marker[iRow] = -1;
2510    mergeLink[iRow] = -1;
2511    link_[iRow] = -1; // not needed but makes debugging easier
2512  }
2513  int start = 0;
2514  int end = 0;
2515  choleskyStart_[0] = 0;
2516
2517  for (iRow = 0; iRow < numberRows_; iRow++) {
2518    int nz = 0;
2519    int merge = mergeLink[iRow];
2520    bool marked = false;
2521    if (merge < 0)
2522      marker[iRow] = iRow;
2523    else
2524      marker[iRow] = merge;
2525    start = end;
2526    int startSub = start;
2527    link_[iRow] = numberRows_;
2528    CoinBigIndex j;
2529    for (j = Astart[iRow]; j < Astart[iRow + 1]; j++) {
2530      int kRow = Arow[j];
2531      int k = iRow;
2532      int linked = link_[iRow];
2533#ifndef NDEBUG
2534      int count = 0;
2535#endif
2536      while (linked <= kRow) {
2537        k = linked;
2538        linked = link_[k];
2539#ifndef NDEBUG
2540        count++;
2541        assert(count < numberRows_);
2542#endif
2543      }
2544      nz++;
2545      link_[k] = kRow;
2546      link_[kRow] = linked;
2547      if (marker[kRow] != marker[iRow])
2548        marked = true;
2549    }
2550    bool reuse = false;
2551    // Check if we can re-use indices
2552    if (!marked && merge >= 0 && mergeLink[merge] < 0) {
2553      // can re-use all
2554      startSub = indexStart_[merge] + 1;
2555      nz = choleskyStart_[merge + 1] - (choleskyStart_[merge] + 1);
2556      reuse = true;
2557    } else {
2558      // See if we can re-use any
2559      int k = mergeLink[iRow];
2560      int maxLength = 0;
2561      while (k >= 0) {
2562        int length = choleskyStart_[k + 1] - (choleskyStart_[k] + 1);
2563        int start = indexStart_[k] + 1;
2564        int stop = start + length;
2565        if (length > maxLength) {
2566          maxLength = length;
2567          startSub = start;
2568        }
2569        int linked = iRow;
2570        for (CoinBigIndex j = start; j < stop; j++) {
2571          int kRow = choleskyRow_[j];
2572          int kk = linked;
2573          linked = link_[kk];
2574          while (linked < kRow) {
2575            kk = linked;
2576            linked = link_[kk];
2577          }
2578          if (linked != kRow) {
2579            nz++;
2580            link_[kk] = kRow;
2581            link_[kRow] = linked;
2582            linked = kRow;
2583          }
2584        }
2585        k = mergeLink[k];
2586      }
2587      if (nz == maxLength)
2588        reuse = true; // can re-use
2589    }
2590    //reuse=false; //temp
2591    if (!reuse) {
2592      end += nz;
2593      startSub = start;
2594      int kRow = iRow;
2595      for (int j = start; j < end; j++) {
2596        kRow = link_[kRow];
2597        choleskyRow_[j] = kRow;
2598        assert(kRow < numberRows_);
2599        marker[kRow] = iRow;
2600      }
2601      marker[iRow] = iRow;
2602    }
2603    indexStart_[iRow] = startSub;
2604    choleskyStart_[iRow + 1] = choleskyStart_[iRow] + nz;
2605    if (nz > 1) {
2606      int kRow = choleskyRow_[startSub];
2607      mergeLink[iRow] = mergeLink[kRow];
2608      mergeLink[kRow] = iRow;
2609    }
2610    // should not be needed
2611    //std::sort(choleskyRow_+indexStart_[iRow]
2612    //      ,choleskyRow_+indexStart_[iRow]+nz);
2613    //#define CLP_DEBUG
2614#ifdef CLP_DEBUG
2615    int last = -1;
2616    for (j = indexStart_[iRow]; j < indexStart_[iRow] + nz; j++) {
2617      int kRow = choleskyRow_[j];
2618      assert(kRow > last);
2619      last = kRow;
2620    }
2621#endif
2622  }
2623  sizeFactor_ = choleskyStart_[numberRows_];
2624  sizeIndex_ = start;
2625  // find dense segment here
2626  int numberleft = numberRows_;
2627  for (iRow = 0; iRow < numberRows_; iRow++) {
2628    CoinBigIndex left = sizeFactor_ - choleskyStart_[iRow];
2629    double n = numberleft;
2630    double threshold = n * (n - 1.0) * 0.5 * goDense_;
2631    if (left >= threshold)
2632      break;
2633    numberleft--;
2634  }
2635  //iRow=numberRows_;
2636  int nDense = numberRows_ - iRow;
2637#define DENSE_THRESHOLD 8
2638  // don't do if dense columns
2639  if (nDense >= DENSE_THRESHOLD && !dense_) {
2640    COIN_DETAIL_PRINT(printf("Going dense for last %d rows\n", nDense));
2641    // make sure we don't disturb any indices
2642    int k = 0;
2643    for (int jRow = 0; jRow < iRow; jRow++) {
2644      int nz = choleskyStart_[jRow + 1] - choleskyStart_[jRow];
2645      k = CoinMax(k, indexStart_[jRow] + nz);
2646    }
2647    indexStart_[iRow] = k;
2648    int j;
2649    for (j = iRow + 1; j < numberRows_; j++) {
2650      choleskyRow_[k++] = j;
2651      indexStart_[j] = k;
2652    }
2653    sizeIndex_ = k;
2654    assert(k <= sizeFactor_); // can't happen with any reasonable defaults
2655    k = choleskyStart_[iRow];
2656    for (j = iRow + 1; j <= numberRows_; j++) {
2657      k += numberRows_ - j;
2658      choleskyStart_[j] = k;
2659    }
2660    // allow for blocked dense
2661    ClpCholeskyDense dense;
2662    sizeFactor_ = choleskyStart_[iRow] + dense.space(nDense);
2663    firstDense_ = iRow;
2664    if (doKKT_) {
2665      // redo permute so negative ones first
2666      int putN = firstDense_;
2667      int putP = 0;
2668      int numberRowsModel = model_->numberRows();
2669      int numberColumns = model_->numberColumns();
2670      int numberTotal = numberColumns + numberRowsModel;
2671      for (iRow = firstDense_; iRow < numberRows_; iRow++) {
2672        int originalRow = permute_[iRow];
2673        if (originalRow < numberTotal)
2674          permute_[putN++] = originalRow;
2675        else
2676          permuteInverse_[putP++] = originalRow;
2677      }
2678      for (iRow = putN; iRow < numberRows_; iRow++) {
2679        permute_[iRow] = permuteInverse_[iRow - putN];
2680      }
2681      for (iRow = 0; iRow < numberRows_; iRow++) {
2682        permuteInverse_[permute_[iRow]] = iRow;
2683      }
2684    }
2685  }
2686  // Clean up clique info
2687  for (iRow = 0; iRow < numberRows_; iRow++)
2688    clique_[iRow] = 0;
2689  int lastClique = -1;
2690  bool inClique = false;
2691  for (iRow = 1; iRow < firstDense_; iRow++) {
2692    int sizeLast = choleskyStart_[iRow] - choleskyStart_[iRow - 1];
2693    int sizeThis = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
2694    if (indexStart_[iRow] == indexStart_[iRow - 1] + 1 && sizeThis == sizeLast - 1 && sizeThis) {
2695      // in clique
2696      if (!inClique) {
2697        inClique = true;
2698        lastClique = iRow - 1;
2699      }
2700    } else if (inClique) {
2701      int sizeClique = iRow - lastClique;
2702      for (int i = lastClique; i < iRow; i++) {
2703        clique_[i] = sizeClique;
2704        sizeClique--;
2705      }
2706      inClique = false;
2707    }
2708  }
2709  if (inClique) {
2710    int sizeClique = iRow - lastClique;
2711    for (int i = lastClique; i < iRow; i++) {
2712      clique_[i] = sizeClique;
2713      sizeClique--;
2714    }
2715  }
2716  //for (iRow=0;iRow<numberRows_;iRow++)
2717  //clique_[iRow]=0;
2718}
2719/* Factorize - filling in rowsDropped and returning number dropped */
2720int ClpCholeskyBase::factorize(const CoinWorkDouble *diagonal, int *rowsDropped)
2721{
2722  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
2723  const int *columnLength = model_->clpMatrix()->getVectorLengths();
2724  const int *row = model_->clpMatrix()->getIndices();
2725  const double *element = model_->clpMatrix()->getElements();
2726  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
2727  const int *rowLength = rowCopy_->getVectorLengths();
2728  const int *column = rowCopy_->getIndices();
2729  const double *elementByRow = rowCopy_->getElements();
2730  int numberColumns = model_->clpMatrix()->getNumCols();
2731  //perturbation
2732  CoinWorkDouble perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
2733  //perturbation=perturbation*perturbation*100000000.0;
2734  if (perturbation > 1.0) {
2735#ifdef COIN_DEVELOP
2736    //if (model_->model()->logLevel()&4)
2737    std::cout << "large perturbation " << perturbation << std::endl;
2738#endif
2739    perturbation = CoinSqrt(perturbation);
2740    perturbation = 1.0;
2741  }
2742  int iRow;
2743  int iColumn;
2744  longDouble *work = workDouble_;
2745  CoinZeroN(work, numberRows_);
2746  int newDropped = 0;
2747  CoinWorkDouble largest = 1.0;
2748  CoinWorkDouble smallest = COIN_DBL_MAX;
2749  int numberDense = 0;
2750  if (!doKKT_) {
2751    const CoinWorkDouble *diagonalSlack = diagonal + numberColumns;
2752    if (dense_)
2753      numberDense = dense_->numberRows();
2754    if (whichDense_) {
2755      longDouble *denseDiagonal = dense_->diagonal();
2756      longDouble *dense = denseColumn_;
2757      int iDense = 0;
2758      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2759        if (whichDense_[iColumn]) {
2760          CoinZeroN(dense, numberRows_);
2761          CoinBigIndex start = columnStart[iColumn];
2762          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2763          if (diagonal[iColumn]) {
2764            denseDiagonal[iDense++] = 1.0 / diagonal[iColumn];
2765            for (CoinBigIndex j = start; j < end; j++) {
2766              int jRow = row[j];
2767              dense[jRow] = element[j];
2768            }
2769          } else {
2770            denseDiagonal[iDense++] = 1.0;
2771          }
2772          dense += numberRows_;
2773        }
2774      }
2775    }
2776    CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
2777    delta2 *= delta2;
2778    // largest in initial matrix
2779    CoinWorkDouble largest2 = 1.0e-20;
2780    for (iRow = 0; iRow < numberRows_; iRow++) {
2781      longDouble *put = sparseFactor_ + choleskyStart_[iRow];
2782      int *which = choleskyRow_ + indexStart_[iRow];
2783      int iOriginalRow = permute_[iRow];
2784      int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
2785      if (!rowLength[iOriginalRow])
2786        rowsDropped_[iOriginalRow] = 1;
2787      if (!rowsDropped_[iOriginalRow]) {
2788        CoinBigIndex startRow = rowStart[iOriginalRow];
2789        CoinBigIndex endRow = rowStart[iOriginalRow] + rowLength[iOriginalRow];
2790        work[iRow] = diagonalSlack[iOriginalRow] + delta2;
2791        for (CoinBigIndex k = startRow; k < endRow; k++) {
2792          int iColumn = column[k];
2793          if (!whichDense_ || !whichDense_[iColumn]) {
2794            CoinBigIndex start = columnStart[iColumn];
2795            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
2796            CoinWorkDouble multiplier = diagonal[iColumn] * elementByRow[k];
2797            for (CoinBigIndex j = start; j < end; j++) {
2798              int jRow = row[j];
2799              int jNewRow = permuteInverse_[jRow];
2800              if (jNewRow >= iRow && !rowsDropped_[jRow]) {
2801                CoinWorkDouble value = element[j] * multiplier;
2802                work[jNewRow] += value;
2803              }
2804            }
2805          }
2806        }
2807        diagonal_[iRow] = work[iRow];
2808        largest2 = CoinMax(largest2, CoinAbs(work[iRow]));
2809        work[iRow] = 0.0;
2810        int j;
2811        for (j = 0; j < number; j++) {
2812          int jRow = which[j];
2813          put[j] = work[jRow];
2814          largest2 = CoinMax(largest2, CoinAbs(work[jRow]));
2815          work[jRow] = 0.0;
2816        }
2817      } else {
2818        // dropped
2819        diagonal_[iRow] = 1.0;
2820        int j;
2821        for (j = 1; j < number; j++) {
2822          put[j] = 0.0;
2823        }
2824      }
2825    }
2826    //check sizes
2827    largest2 *= 1.0e-20;
2828    largest = CoinMin(largest2, CHOL_SMALL_VALUE);
2829    int numberDroppedBefore = 0;
2830    for (iRow = 0; iRow < numberRows_; iRow++) {
2831      int dropped = rowsDropped_[iRow];
2832      // Move to int array
2833      rowsDropped[iRow] = dropped;
2834      if (!dropped) {
2835        CoinWorkDouble diagonal = diagonal_[iRow];
2836        if (diagonal > largest2) {
2837          diagonal_[iRow] = diagonal + perturbation;
2838        } else {
2839          diagonal_[iRow] = diagonal + perturbation;
2840          rowsDropped[iRow] = 2;
2841          numberDroppedBefore++;
2842          //printf("dropped - small diagonal %g\n",diagonal);
2843        }
2844      }
2845    }
2846    doubleParameters_[10] = CoinMax(1.0e-20, largest);
2847    integerParameters_[20] = 0;
2848    doubleParameters_[3] = 0.0;
2849    doubleParameters_[4] = COIN_DBL_MAX;
2850    integerParameters_[34] = 0; // say all must be positive
2851    factorizePart2(rowsDropped);
2852    newDropped = integerParameters_[20] + numberDroppedBefore;
2853    largest = doubleParameters_[3];
2854    smallest = doubleParameters_[4];
2855    if (model_->messageHandler()->logLevel() > 1)
2856      std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
2857    choleskyCondition_ = largest / smallest;
2858    if (whichDense_) {
2859      int i;
2860      for (i = 0; i < numberRows_; i++) {
2861        assert(diagonal_[i] >= 0.0);
2862        diagonal_[i] = CoinSqrt(diagonal_[i]);
2863      }
2864      // Update dense columns (just L)
2865      // Zero out dropped rows
2866      for (i = 0; i < numberDense; i++) {
2867        longDouble *a = denseColumn_ + i * numberRows_;
2868        for (int j = 0; j < numberRows_; j++) {
2869          if (rowsDropped[j])
2870            a[j] = 0.0;
2871        }
2872        for (i = 0; i < numberRows_; i++) {
2873          int iRow = permute_[i];
2874          workDouble_[i] = a[iRow];
2875        }
2876        for (i = 0; i < numberRows_; i++) {
2877          CoinWorkDouble value = workDouble_[i];
2878          CoinBigIndex offset = indexStart_[i] - choleskyStart_[i];
2879          CoinBigIndex j;
2880          for (j = choleskyStart_[i]; j < choleskyStart_[i + 1]; j++) {
2881            int iRow = choleskyRow_[j + offset];
2882            workDouble_[iRow] -= sparseFactor_[j] * value;
2883          }
2884        }
2885        for (i = 0; i < numberRows_; i++) {
2886          int iRow = permute_[i];
2887          a[iRow] = workDouble_[i] * diagonal_[i];
2888        }
2889      }
2890      dense_->resetRowsDropped();
2891      longDouble *denseBlob = dense_->aMatrix();
2892      longDouble *denseDiagonal = dense_->diagonal();
2893      // Update dense matrix
2894      for (i = 0; i < numberDense; i++) {
2895        const longDouble *a = denseColumn_ + i * numberRows_;
2896        // do diagonal
2897        CoinWorkDouble value = denseDiagonal[i];
2898        const longDouble *b = denseColumn_ + i * numberRows_;
2899        for (int k = 0; k < numberRows_; k++)
2900          value += a[k] * b[k];
2901        denseDiagonal[i] = value;
2902        for (int j = i + 1; j < numberDense; j++) {
2903          CoinWorkDouble value = 0.0;
2904          const longDouble *b = denseColumn_ + j * numberRows_;
2905          for (int k = 0; k < numberRows_; k++)
2906            value += a[k] * b[k];
2907          *denseBlob = value;
2908          denseBlob++;
2909        }
2910      }
2911      // dense cholesky (? long double)
2912      int *dropped = new int[numberDense];
2913      dense_->factorizePart2(dropped);
2914      delete[] dropped;
2915    }
2916    // try allowing all every time
2917    //printf("trying ?\n");
2918    //for (iRow=0;iRow<numberRows_;iRow++) {
2919    //rowsDropped[iRow]=0;
2920    //rowsDropped_[iRow]=0;
2921    //}
2922    bool cleanCholesky;
2923    //if (model_->numberIterations()<20||(model_->numberIterations()&1)==0)
2924    if (model_->numberIterations() < 2000)
2925      cleanCholesky = true;
2926    else
2927      cleanCholesky = false;
2928    if (cleanCholesky) {
2929      //drop fresh makes some formADAT easier
2930      if (newDropped || numberRowsDropped_) {
2931        newDropped = 0;
2932        for (int i = 0; i < numberRows_; i++) {
2933          char dropped = static_cast< char >(rowsDropped[i]);
2934          rowsDropped_[i] = dropped;
2935          rowsDropped_[i] = 0;
2936          if (dropped == 2) {
2937            //dropped this time
2938            rowsDropped[newDropped++] = i;
2939            rowsDropped_[i] = 0;
2940          }
2941        }
2942        numberRowsDropped_ = newDropped;
2943        newDropped = -(2 + newDropped);
2944      }
2945    } else {
2946      if (newDropped) {
2947        newDropped = 0;
2948        for (int i = 0; i < numberRows_; i++) {
2949          char dropped = static_cast< char >(rowsDropped[i]);
2950          rowsDropped_[i] = dropped;
2951          if (dropped == 2) {
2952            //dropped this time
2953            rowsDropped[newDropped++] = i;
2954            rowsDropped_[i] = 1;
2955          }
2956        }
2957      }
2958      numberRowsDropped_ += newDropped;
2959      if (numberRowsDropped_ && 0) {
2960        std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)";
2961        if (newDropped) {
2962          std::cout << " ( " << newDropped << " dropped this time)";
2963        }
2964        std::cout << std::endl;
2965      }
2966    }
2967  } else {
2968    //KKT
2969    CoinPackedMatrix *quadratic = NULL;
2970    ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(model_->objectiveAsObject()));
2971    if (quadraticObj)
2972      quadratic = quadraticObj->quadraticObjective();
2973    // matrix
2974    int numberRowsModel = model_->numberRows();
2975    int numberColumns = model_->numberColumns();
2976    int numberTotal = numberColumns + numberRowsModel;
2977    // temp
2978    bool permuted = false;
2979    for (iRow = 0; iRow < numberRows_; iRow++) {
2980      if (permute_[iRow] != iRow) {
2981        permuted = true;
2982        break;
2983      }
2984    }
2985    // but fake it
2986    for (iRow = 0; iRow < numberRows_; iRow++) {
2987      //permute_[iRow]=iRow; // force no permute
2988      //permuteInverse_[permute_[iRow]]=iRow;
2989    }
2990    if (permuted) {
2991      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
2992      delta2 *= delta2;
2993      // Need to permute - ugly
2994      if (!quadratic) {
2995        for (iRow = 0; iRow < numberRows_; iRow++) {
2996          longDouble *put = sparseFactor_ + choleskyStart_[iRow];
2997          int *which = choleskyRow_ + indexStart_[iRow];
2998          int iOriginalRow = permute_[iRow];
2999          if (iOriginalRow < numberColumns) {
3000            iColumn = iOriginalRow;
3001            CoinWorkDouble value = diagonal[iColumn];
3002            if (CoinAbs(value) > 1.0e-100) {
3003              value = 1.0 / value;
3004              largest = CoinMax(largest, CoinAbs(value));
3005              diagonal_[iRow] = -value;
3006              CoinBigIndex start = columnStart[iColumn];
3007              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
3008              for (CoinBigIndex j = start; j < end; j++) {
3009                int kRow = row[j] + numberTotal;
3010                kRow = permuteInverse_[kRow];
3011                if (kRow > iRow) {
3012                  work[kRow] = element[j];
3013                  largest = CoinMax(largest, CoinAbs(element[j]));
3014                }
3015              }
3016            } else {
3017              diagonal_[iRow] = -value;
3018            }
3019          } else if (iOriginalRow < numberTotal) {
3020            CoinWorkDouble value = diagonal[iOriginalRow];
3021            if (CoinAbs(value) > 1.0e-100) {
3022              value = 1.0 / value;
3023              largest = CoinMax(largest, CoinAbs(value));
3024            } else {
3025              value = 1.0e100;
3026            }
3027            diagonal_[iRow] = -value;
3028            int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
3029            if (kRow > iRow)
3030              work[kRow] = -1.0;
3031          } else {
3032            diagonal_[iRow] = delta2;
3033            int kRow = iOriginalRow - numberTotal;
3034            CoinBigIndex start = rowStart[kRow];
3035            CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
3036            for (CoinBigIndex j = start; j < end; j++) {
3037              int jRow = column[j];
3038              int jNewRow = permuteInverse_[jRow];
3039              if (jNewRow > iRow) {
3040                work[jNewRow] = elementByRow[j];
3041                largest = CoinMax(largest, CoinAbs(elementByRow[j]));
3042              }
3043            }
3044            // slack - should it be permute
3045            kRow = permuteInverse_[kRow + numberColumns];
3046            if (kRow > iRow)
3047              work[kRow] = -1.0;
3048          }
3049          CoinBigIndex j;
3050          int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
3051          for (j = 0; j < number; j++) {
3052            int jRow = which[j];
3053            put[j] = work[jRow];
3054            work[jRow] = 0.0;
3055          }
3056        }
3057      } else {
3058        // quadratic
3059        const int *columnQuadratic = quadratic->getIndices();
3060        const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
3061        const int *columnQuadraticLength = quadratic->getVectorLengths();
3062        const double *quadraticElement = quadratic->getElements();
3063        for (iRow = 0; iRow < numberRows_; iRow++) {
3064          longDouble *put = sparseFactor_ + choleskyStart_[iRow];
3065          int *which = choleskyRow_ + indexStart_[iRow];
3066          int iOriginalRow = permute_[iRow];
3067          if (iOriginalRow < numberColumns) {
3068            CoinBigIndex j;
3069            iColumn = iOriginalRow;
3070            CoinWorkDouble value = diagonal[iColumn];
3071            if (CoinAbs(value) > 1.0e-100) {
3072              value = 1.0 / value;
3073              for (j = columnQuadraticStart[iColumn];
3074                   j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
3075                int jColumn = columnQuadratic[j];
3076                int jNewColumn = permuteInverse_[jColumn];
3077                if (jNewColumn > iRow) {
3078                  work[jNewColumn] = -quadraticElement[j];
3079                } else if (iColumn == jColumn) {
3080                  value += quadraticElement[j];
3081                }
3082              }
3083              largest = CoinMax(largest, CoinAbs(value));
3084              diagonal_[iRow] = -value;
3085              CoinBigIndex start = columnStart[iColumn];
3086              CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
3087              for (j = start; j < end; j++) {
3088                int kRow = row[j] + numberTotal;
3089                kRow = permuteInverse_[kRow];
3090                if (kRow > iRow) {
3091                  work[kRow] = element[j];
3092                  largest = CoinMax(largest, CoinAbs(element[j]));
3093                }
3094              }
3095            } else {
3096              diagonal_[iRow] = -value;
3097            }
3098          } else if (iOriginalRow < numberTotal) {
3099            CoinWorkDouble value = diagonal[iOriginalRow];
3100            if (CoinAbs(value) > 1.0e-100) {
3101              value = 1.0 / value;
3102              largest = CoinMax(largest, CoinAbs(value));
3103            } else {
3104              value = 1.0e100;
3105            }
3106            diagonal_[iRow] = -value;
3107            int kRow = permuteInverse_[iOriginalRow + numberRowsModel];
3108            if (kRow > iRow)
3109              work[kRow] = -1.0;
3110          } else {
3111            diagonal_[iRow] = delta2;
3112            int kRow = iOriginalRow - numberTotal;
3113            CoinBigIndex start = rowStart[kRow];
3114            CoinBigIndex end = rowStart[kRow] + rowLength[kRow];
3115            for (CoinBigIndex j = start; j < end; j++) {
3116              int jRow = column[j];
3117              int jNewRow = permuteInverse_[jRow];
3118              if (jNewRow > iRow) {
3119                work[jNewRow] = elementByRow[j];
3120                largest = CoinMax(largest, CoinAbs(elementByRow[j]));
3121              }
3122            }
3123            // slack - should it be permute
3124            kRow = permuteInverse_[kRow + numberColumns];
3125            if (kRow > iRow)
3126              work[kRow] = -1.0;
3127          }
3128          CoinBigIndex j;
3129          int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
3130          for (j = 0; j < number; j++) {
3131            int jRow = which[j];
3132            put[j] = work[jRow];
3133            work[jRow] = 0.0;
3134          }
3135          for (j = 0; j < numberRows_; j++)
3136            assert(!work[j]);
3137        }
3138      }
3139    } else {
3140      if (!quadratic) {
3141        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3142          longDouble *put = sparseFactor_ + choleskyStart_[iColumn];
3143          int *which = choleskyRow_ + indexStart_[iColumn];
3144          CoinWorkDouble value = diagonal[iColumn];
3145          if (CoinAbs(value) > 1.0e-100) {
3146            value = 1.0 / value;
3147            largest = CoinMax(largest, CoinAbs(value));
3148            diagonal_[iColumn] = -value;
3149            CoinBigIndex start = columnStart[iColumn];
3150            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
3151            for (CoinBigIndex j = start; j < end; j++) {
3152              //choleskyRow_[numberElements]=row[j]+numberTotal;
3153              //sparseFactor_[numberElements++]=element[j];
3154              work[row[j] + numberTotal] = element[j];
3155              largest = CoinMax(largest, CoinAbs(element[j]));
3156            }
3157          } else {
3158            diagonal_[iColumn] = -value;
3159          }
3160          CoinBigIndex j;
3161          int number = choleskyStart_[iColumn + 1] - choleskyStart_[iColumn];
3162          for (j = 0; j < number; j++) {
3163            int jRow = which[j];
3164            put[j] = work[jRow];
3165            work[jRow] = 0.0;
3166          }
3167        }
3168      } else {
3169        // Quadratic
3170        const int *columnQuadratic = quadratic->getIndices();
3171        const CoinBigIndex *columnQuadraticStart = quadratic->getVectorStarts();
3172        const int *columnQuadraticLength = quadratic->getVectorLengths();
3173        const double *quadraticElement = quadratic->getElements();
3174        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3175          longDouble *put = sparseFactor_ + choleskyStart_[iColumn];
3176          int *which = choleskyRow_ + indexStart_[iColumn];
3177          int number = choleskyStart_[iColumn + 1] - choleskyStart_[iColumn];
3178          CoinWorkDouble value = diagonal[iColumn];
3179          CoinBigIndex j;
3180          if (CoinAbs(value) > 1.0e-100) {
3181            value = 1.0 / value;
3182            for (j = columnQuadraticStart[iColumn];
3183                 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
3184              int jColumn = columnQuadratic[j];
3185              if (jColumn > iColumn) {
3186                work[jColumn] = -quadraticElement[j];
3187              } else if (iColumn == jColumn) {
3188                value += quadraticElement[j];
3189              }
3190            }
3191            largest = CoinMax(largest, CoinAbs(value));
3192            diagonal_[iColumn] = -value;
3193            CoinBigIndex start = columnStart[iColumn];
3194            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
3195            for (j = start; j < end; j++) {
3196              work[row[j] + numberTotal] = element[j];
3197              largest = CoinMax(largest, CoinAbs(element[j]));
3198            }
3199            for (j = 0; j < number; j++) {
3200              int jRow = which[j];
3201              put[j] = work[jRow];
3202              work[jRow] = 0.0;
3203            }
3204          } else {
3205            value = 1.0e100;
3206            diagonal_[iColumn] = -value;
3207            for (j = 0; j < number; j++) {
3208              int jRow = which[j];
3209              put[j] = work[jRow];
3210            }
3211          }
3212        }
3213      }
3214      // slacks
3215      for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) {
3216        longDouble *put = sparseFactor_ + choleskyStart_[iColumn];
3217        int *which = choleskyRow_ + indexStart_[iColumn];
3218        CoinWorkDouble value = diagonal[iColumn];
3219        if (CoinAbs(value) > 1.0e-100) {
3220          value = 1.0 / value;
3221          largest = CoinMax(largest, CoinAbs(value));
3222        } else {
3223          value = 1.0e100;
3224        }
3225        diagonal_[iColumn] = -value;
3226        work[iColumn - numberColumns + numberTotal] = -1.0;
3227        CoinBigIndex j;
3228        int number = choleskyStart_[iColumn + 1] - choleskyStart_[iColumn];
3229        for (j = 0; j < number; j++) {
3230          int jRow = which[j];
3231          put[j] = work[jRow];
3232          work[jRow] = 0.0;
3233        }
3234      }
3235      // Finish diagonal
3236      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
3237      delta2 *= delta2;
3238      for (iRow = 0; iRow < numberRowsModel; iRow++) {
3239        longDouble *put = sparseFactor_ + choleskyStart_[iRow + numberTotal];
3240        diagonal_[iRow + numberTotal] = delta2;
3241        CoinBigIndex j;
3242        int number = choleskyStart_[iRow + numberTotal + 1] - choleskyStart_[iRow + numberTotal];
3243        for (j = 0; j < number; j++) {
3244          put[j] = 0.0;
3245        }
3246      }
3247    }
3248    //check sizes
3249    largest *= 1.0e-20;
3250    largest = CoinMin(largest, CHOL_SMALL_VALUE);
3251    doubleParameters_[10] = CoinMax(1.0e-20, largest);
3252    integerParameters_[20] = 0;
3253    doubleParameters_[3] = 0.0;
3254    doubleParameters_[4] = COIN_DBL_MAX;
3255    // Set up LDL cutoff
3256    integerParameters_[34] = numberTotal;
3257    // KKT
3258    int *rowsDropped2 = new int[numberRows_];
3259    CoinZeroN(rowsDropped2, numberRows_);
3260    factorizePart2(rowsDropped2);
3261    newDropped = integerParameters_[20];
3262    largest = doubleParameters_[3];
3263    smallest = doubleParameters_[4];
3264    if (model_->messageHandler()->logLevel() > 1)
3265      std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
3266    choleskyCondition_ = largest / smallest;
3267    // Should save adjustments in ..R_
3268    int n1 = 0, n2 = 0;
3269    CoinWorkDouble *primalR = model_->primalR();
3270    CoinWorkDouble *dualR = model_->dualR();
3271    for (iRow = 0; iRow < numberTotal; iRow++) {
3272      if (rowsDropped2[iRow]) {
3273        n1++;
3274        //printf("row region1 %d dropped\n",iRow);
3275        //rowsDropped_[iRow]=1;
3276        rowsDropped_[iRow] = 0;
3277        primalR[iRow] = doubleParameters_[20];
3278      } else {
3279        rowsDropped_[iRow] = 0;
3280        primalR[iRow] = 0.0;
3281      }
3282    }
3283    for (; iRow < numberRows_; iRow++) {
3284      if (rowsDropped2[iRow]) {
3285        n2++;
3286        //printf("row region2 %d dropped\n",iRow);
3287        //rowsDropped_[iRow]=1;
3288        rowsDropped_[iRow] = 0;
3289        dualR[iRow - numberTotal] = doubleParameters_[34];
3290      } else {
3291        rowsDropped_[iRow] = 0;
3292        dualR[iRow - numberTotal] = 0.0;
3293      }
3294    }
3295    delete[] rowsDropped2;
3296  }
3297  status_ = 0;
3298  return newDropped;
3299}
3300/* Factorize - filling in rowsDropped and returning number dropped
3301   in integerParam.
3302*/
3303void ClpCholeskyBase::factorizePart2(int *rowsDropped)
3304{
3305  CoinWorkDouble largest = doubleParameters_[3];
3306  CoinWorkDouble smallest = doubleParameters_[4];
3307  // probably done before
3308  largest = 0.0;
3309  smallest = COIN_DBL_MAX;
3310  double dropValue = doubleParameters_[10];
3311  int firstPositive = integerParameters_[34];
3312  longDouble *d = ClpCopyOfArray(diagonal_, numberRows_);
3313  int iRow;
3314  // minimum size before clique done
3315  //#define MINCLIQUE INT_MAX
3316#define MINCLIQUE 3
3317  longDouble *work = workDouble_;
3318  int *first = workInteger_;
3319
3320  for (iRow = 0; iRow < numberRows_; iRow++) {
3321    link_[iRow] = -1;
3322    work[iRow] = 0.0;
3323    first[iRow] = choleskyStart_[iRow];
3324  }
3325
3326  int lastClique = -1;
3327  bool inClique = false;
3328  bool newClique = false;
3329  bool endClique = false;
3330  int lastRow = 0;
3331  int nextRow2 = -1;
3332
3333  for (iRow = 0; iRow < firstDense_ + 1; iRow++) {
3334    if (iRow < firstDense_) {
3335      endClique = false;
3336      if (clique_[iRow] > 0) {
3337        // this is in a clique
3338        inClique = true;
3339        if (clique_[iRow] > lastClique) {
3340          // new Clique
3341          newClique = true;
3342          // If we have clique going then signal to do old one
3343          endClique = (lastClique > 0);
3344        } else {
3345          // Still in clique
3346          newClique = false;
3347        }
3348      } else {
3349        // not in clique
3350        inClique = false;
3351        newClique = false;
3352        // If we have clique going then signal to do old one
3353        endClique = (lastClique > 0);
3354      }
3355      lastClique = clique_[iRow];
3356    } else if (inClique) {
3357      // Finish off
3358      endClique = true;
3359    } else {
3360      break;
3361    }
3362    if (endClique) {
3363      // We have just finished updating a clique - do block pivot and clean up
3364      int jRow;
3365      for (jRow = lastRow; jRow < iRow; jRow++) {
3366        int jCount = jRow - lastRow;
3367        CoinWorkDouble diagonalValue = diagonal_[jRow];
3368        CoinBigIndex start = choleskyStart_[jRow];
3369        CoinBigIndex end = choleskyStart_[jRow + 1];
3370        for (int kRow = lastRow; kRow < jRow; kRow++) {
3371          jCount--;
3372          CoinBigIndex get = choleskyStart_[kRow] + jCount;
3373          CoinWorkDouble a_jk = sparseFactor_[get];
3374          CoinWorkDouble value1 = d[kRow] * a_jk;
3375          diagonalValue -= a_jk * value1;
3376          for (CoinBigIndex j = start; j < end; j++)
3377            sparseFactor_[j] -= value1 * sparseFactor_[++get];
3378        }
3379        // check
3380        int originalRow = permute_[jRow];
3381        if (originalRow < firstPositive) {
3382          // must be negative
3383          if (diagonalValue <= -dropValue) {
3384            smallest = CoinMin(smallest, -diagonalValue);
3385            largest = CoinMax(largest, -diagonalValue);
3386            d[jRow] = diagonalValue;
3387            diagonalValue = 1.0 / diagonalValue;
3388          } else {
3389            rowsDropped[originalRow] = 2;
3390            d[jRow] = -1.0e100;
3391            diagonalValue = 0.0;
3392            integerParameters_[20]++;
3393          }
3394        } else {
3395          // must be positive
3396          if (diagonalValue >= dropValue) {
3397            smallest = CoinMin(smallest, diagonalValue);
3398            largest = CoinMax(largest, diagonalValue);
3399            d[jRow] = diagonalValue;
3400            diagonalValue = 1.0 / diagonalValue;
3401          } else {
3402            rowsDropped[originalRow] = 2;
3403            d[jRow] = 1.0e100;
3404            diagonalValue = 0.0;
3405            integerParameters_[20]++;
3406          }
3407        }
3408        diagonal_[jRow] = diagonalValue;
3409        for (CoinBigIndex j = start; j < end; j++) {
3410          sparseFactor_[j] *= diagonalValue;
3411        }
3412      }
3413      if (nextRow2 >= 0) {
3414        for (jRow = lastRow; jRow < iRow - 1; jRow++) {
3415          link_[jRow] = jRow + 1;
3416        }
3417        link_[iRow - 1] = link_[nextRow2];
3418        link_[nextRow2] = lastRow;
3419      }
3420    }
3421    if (iRow == firstDense_)
3422      break; // we were just cleaning up
3423    if (newClique) {
3424      // initialize new clique
3425      lastRow = iRow;
3426    }
3427    // for each column L[*,kRow] that affects L[*,iRow]
3428    CoinWorkDouble diagonalValue = diagonal_[iRow];
3429    int nextRow = link_[iRow];
3430    int kRow = 0;
3431    while (1) {
3432      kRow = nextRow;
3433      if (kRow < 0)
3434        break; // out of loop
3435      nextRow = link_[kRow];
3436      // Modify by outer product of L[*,irow] by L[*,krow] from first
3437      int k = first[kRow];
3438      int end = choleskyStart_[kRow + 1];
3439      assert(k < end);
3440      CoinWorkDouble a_ik = sparseFactor_[k++];
3441      CoinWorkDouble value1 = d[kRow] * a_ik;
3442      // update first
3443      first[kRow] = k;
3444      diagonalValue -= value1 * a_ik;
3445      int offset = indexStart_[kRow] - choleskyStart_[kRow];
3446      if (k < end) {
3447        int jRow = choleskyRow_[k + offset];
3448        if (clique_[kRow] < MINCLIQUE) {
3449          link_[kRow] = link_[jRow];
3450          link_[jRow] = kRow;
3451          for (; k < end; k++) {
3452            int jRow = choleskyRow_[k + offset];
3453            work[jRow] += sparseFactor_[k] * value1;
3454          }
3455        } else {
3456          // Clique
3457          CoinBigIndex currentIndex = k + offset;
3458          int linkSave = link_[jRow];
3459          link_[jRow] = kRow;
3460          work[kRow] = value1; // ? or a_jk
3461          int last = kRow + clique_[kRow];
3462          for (int kkRow = kRow + 1; kkRow < last; kkRow++) {
3463            CoinBigIndex j = first[kkRow];
3464            //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
3465            CoinWorkDouble a = sparseFactor_[j];
3466            CoinWorkDouble dValue = d[kkRow] * a;
3467            diagonalValue -= a * dValue;
3468            work[kkRow] = dValue;
3469            first[kkRow]++;
3470            link_[kkRow - 1] = kkRow;
3471          }
3472          nextRow = link_[last - 1];
3473          link_[last - 1] = linkSave;
3474          int length = static_cast< int >(end - k);
3475          for (int i = 0; i < length; i++) {
3476            int lRow = choleskyRow_[currentIndex++];
3477            CoinWorkDouble t0 = work[lRow];
3478            for (int kkRow = kRow; kkRow < last; kkRow++) {
3479              CoinBigIndex j = first[kkRow] + i;
3480              t0 += work[kkRow] * sparseFactor_[j];
3481            }
3482            work[lRow] = t0;
3483          }
3484        }
3485      }
3486    }
3487    // Now apply
3488    if (inClique) {
3489      // in clique
3490      diagonal_[iRow] = diagonalValue;
3491      int start = choleskyStart_[iRow];
3492      int end = choleskyStart_[iRow + 1];
3493      int currentIndex = indexStart_[iRow];
3494      nextRow2 = -1;
3495      int get = start + clique_[iRow] - 1;
3496      if (get < end) {
3497        nextRow2 = choleskyRow_[currentIndex + get - start];
3498        first[iRow] = get;
3499      }
3500      for (int j = start; j < end; j++) {
3501        int kRow = choleskyRow_[currentIndex++];
3502        sparseFactor_[j] -= work[kRow]; // times?
3503        work[kRow] = 0.0;
3504      }
3505    } else {
3506      // not in clique
3507      int originalRow = permute_[iRow];
3508      if (originalRow < firstPositive) {
3509        // must be negative
3510        if (diagonalValue <= -dropValue) {
3511          smallest = CoinMin(smallest, -diagonalValue);
3512          largest = CoinMax(largest, -diagonalValue);
3513          d[iRow] = diagonalValue;
3514          diagonalValue = 1.0 / diagonalValue;
3515        } else {
3516          rowsDropped[originalRow] = 2;
3517          d[iRow] = -1.0e100;
3518          diagonalValue = 0.0;
3519          integerParameters_[20]++;
3520        }
3521      } else {
3522        // must be positive
3523        if (diagonalValue >= dropValue) {
3524          smallest = CoinMin(smallest, diagonalValue);
3525          largest = CoinMax(largest, diagonalValue);
3526          d[iRow] = diagonalValue;
3527          diagonalValue = 1.0 / diagonalValue;
3528        } else {
3529          rowsDropped[originalRow] = 2;
3530          d[iRow] = 1.0e100;
3531          diagonalValue = 0.0;
3532          integerParameters_[20]++;
3533        }
3534      }
3535      diagonal_[iRow] = diagonalValue;
3536      int offset = indexStart_[iRow] - choleskyStart_[iRow];
3537      int start = choleskyStart_[iRow];
3538      int end = choleskyStart_[iRow + 1];
3539      assert(first[iRow] == start);
3540      if (start < end) {
3541        int nextRow = choleskyRow_[start + offset];
3542        link_[iRow] = link_[nextRow];
3543        link_[nextRow] = iRow;
3544        for (int j = start; j < end; j++) {
3545          int jRow = choleskyRow_[j + offset];
3546          CoinWorkDouble value = sparseFactor_[j] - work[jRow];
3547          work[jRow] = 0.0;
3548          sparseFactor_[j] = diagonalValue * value;
3549        }
3550      }
3551    }
3552  }
3553  if (firstDense_ < numberRows_) {
3554    // do dense
3555    // update dense part
3556    updateDense(d, /*work,*/ first);
3557    ClpCholeskyDense dense;
3558    // just borrow space
3559    int nDense = numberRows_ - firstDense_;
3560    if (doKKT_) {
3561      for (iRow = firstDense_; iRow < numberRows_; iRow++) {
3562        int originalRow = permute_[iRow];
3563        if (originalRow >= firstPositive) {
3564          firstPositive = iRow - firstDense_;
3565          break;
3566        }
3567      }
3568    }
3569    dense.reserveSpace(this, nDense);
3570    int *dropped = new int[nDense];
3571    memset(dropped, 0, nDense * sizeof(int));
3572    dense.setDoubleParameter(3, largest);
3573    dense.setDoubleParameter(4, smallest);
3574    dense.setDoubleParameter(10, dropValue);
3575    dense.setIntegerParameter(20, 0);
3576    dense.setIntegerParameter(34, firstPositive);
3577    dense.setModel(model_);
3578    dense.factorizePart2(dropped);
3579    largest = dense.getDoubleParameter(3);
3580    smallest = dense.getDoubleParameter(4);
3581    integerParameters_[20] += dense.getIntegerParameter(20);
3582    for (iRow = firstDense_; iRow < numberRows_; iRow++) {
3583      int originalRow = permute_[iRow];
3584      rowsDropped[originalRow] = dropped[iRow - firstDense_];
3585    }
3586    delete[] dropped;
3587  }
3588  delete[] d;
3589  doubleParameters_[3] = largest;
3590  doubleParameters_[4] = smallest;
3591  return;
3592}
3593// Updates dense part (broken out for profiling)
3594void ClpCholeskyBase::updateDense(longDouble *d, /*longDouble * work,*/ int *first)
3595{
3596  for (int iRow = 0; iRow < firstDense_; iRow++) {
3597    int start = first[iRow];
3598    int end = choleskyStart_[iRow + 1];
3599    if (start < end) {
3600      int offset = indexStart_[iRow] - choleskyStart_[iRow];
3601      if (clique_[iRow] < 2) {
3602        CoinWorkDouble dValue = d[iRow];
3603        for (int k = start; k < end; k++) {
3604          int kRow = choleskyRow_[k + offset];
3605          assert(kRow >= firstDense_);
3606          CoinWorkDouble a_ik = sparseFactor_[k];
3607          CoinWorkDouble value1 = dValue * a_ik;
3608          diagonal_[kRow] -= value1 * a_ik;
3609          int base = choleskyStart_[kRow] - kRow - 1;
3610          for (int j = k + 1; j < end; j++) {
3611            int jRow = choleskyRow_[j + offset];
3612            CoinWorkDouble a_jk = sparseFactor_[j];
3613            sparseFactor_[base + jRow] -= a_jk * value1;
3614          }
3615        }
3616      } else if (clique_[iRow] < 3) {
3617        // do as pair
3618        CoinWorkDouble dValue0 = d[iRow];
3619        CoinWorkDouble dValue1 = d[iRow + 1];
3620        int offset1 = first[iRow + 1] - start;
3621        // skip row
3622        iRow++;
3623        for (int k = start; k < end; k++) {
3624          int kRow = choleskyRow_[k + offset];
3625          assert(kRow >= firstDense_);
3626          CoinWorkDouble a_ik0 = sparseFactor_[k];
3627          CoinWorkDouble value0 = dValue0 * a_ik0;
3628          CoinWorkDouble a_ik1 = sparseFactor_[k + offset1];
3629          CoinWorkDouble value1 = dValue1 * a_ik1;
3630          diagonal_[kRow] -= value0 * a_ik0 + value1 * a_ik1;
3631          int base = choleskyStart_[kRow] - kRow - 1;
3632          for (int j = k + 1; j < end; j++) {
3633            int jRow = choleskyRow_[j + offset];
3634            CoinWorkDouble a_jk0 = sparseFactor_[j];
3635            CoinWorkDouble a_jk1 = sparseFactor_[j + offset1];
3636            sparseFactor_[base + jRow] -= a_jk0 * value0 + a_jk1 * value1;
3637          }
3638        }
3639#define MANY_REGISTERS
3640#ifdef MANY_REGISTERS
3641      } else if (clique_[iRow] == 3) {
3642#else
3643      } else {
3644#endif
3645        // do as clique
3646        // maybe later get fancy on big cliques and do transpose copy
3647        // seems only worth going to 3 on Intel
3648        CoinWorkDouble dValue0 = d[iRow];
3649        CoinWorkDouble dValue1 = d[iRow + 1];
3650        CoinWorkDouble dValue2 = d[iRow + 2];
3651        // get offsets and skip rows
3652        int offset1 = first[++iRow] - start;
3653        int offset2 = first[++iRow] - start;
3654        for (int k = start; k < end; k++) {
3655          int kRow = choleskyRow_[k + offset];
3656          assert(kRow >= firstDense_);
3657          CoinWorkDouble diagonalValue = diagonal_[kRow];
3658          CoinWorkDouble a_ik0 = sparseFactor_[k];
3659          CoinWorkDouble value0 = dValue0 * a_ik0;
3660          CoinWorkDouble a_ik1 = sparseFactor_[k + offset1];
3661          CoinWorkDouble value1 = dValue1 * a_ik1;
3662          CoinWorkDouble a_ik2 = sparseFactor_[k + offset2];
3663          CoinWorkDouble value2 = dValue2 * a_ik2;
3664          int base = choleskyStart_[kRow] - kRow - 1;
3665          diagonal_[kRow] = diagonalValue - value0 * a_ik0 - value1 * a_ik1 - value2 * a_ik2;
3666          for (int j = k + 1; j < end; j++) {
3667            int jRow = choleskyRow_[j + offset];
3668            CoinWorkDouble a_jk0 = sparseFactor_[j];
3669            CoinWorkDouble a_jk1 = sparseFactor_[j + offset1];
3670            CoinWorkDouble a_jk2 = sparseFactor_[j + offset2];
3671            sparseFactor_[base + jRow] -= a_jk0 * value0 + a_jk1 * value1 + a_jk2 * value2;
3672          }
3673        }
3674#ifdef MANY_REGISTERS
3675      } else {
3676        // do as clique
3677        // maybe later get fancy on big cliques and do transpose copy
3678        // maybe only worth going to 3 on Intel (but may have hidden registers)
3679        CoinWorkDouble dValue0 = d[iRow];
3680        CoinWorkDouble dValue1 = d[iRow + 1];
3681        CoinWorkDouble dValue2 = d[iRow + 2];
3682        CoinWorkDouble dValue3 = d[iRow + 3];
3683        // get offsets and skip rows
3684        int offset1 = first[++iRow] - start;
3685        int offset2 = first[++iRow] - start;
3686        int offset3 = first[++iRow] - start;
3687        for (int k = start; k < end; k++) {
3688          int kRow = choleskyRow_[k + offset];
3689          assert(kRow >= firstDense_);
3690          CoinWorkDouble diagonalValue = diagonal_[kRow];
3691          CoinWorkDouble a_ik0 = sparseFactor_[k];
3692          CoinWorkDouble value0 = dValue0 * a_ik0;
3693          CoinWorkDouble a_ik1 = sparseFactor_[k + offset1];
3694          CoinWorkDouble value1 = dValue1 * a_ik1;
3695          CoinWorkDouble a_ik2 = sparseFactor_[k + offset2];
3696          CoinWorkDouble value2 = dValue2 * a_ik2;
3697          CoinWorkDouble a_ik3 = sparseFactor_[k + offset3];
3698          CoinWorkDouble value3 = dValue3 * a_ik3;
3699          int base = choleskyStart_[kRow] - kRow - 1;
3700          diagonal_[kRow] = diagonalValue - (value0 * a_ik0 + value1 * a_ik1 + value2 * a_ik2 + value3 * a_ik3);
3701          for (int j = k + 1; j < end; j++) {
3702            int jRow = choleskyRow_[j + offset];
3703            CoinWorkDouble a_jk0 = sparseFactor_[j];
3704            CoinWorkDouble a_jk1 = sparseFactor_[j + offset1];
3705            CoinWorkDouble a_jk2 = sparseFactor_[j + offset2];
3706            CoinWorkDouble a_jk3 = sparseFactor_[j + offset3];
3707            sparseFactor_[base + jRow] -= a_jk0 * value0 + a_jk1 * value1 + a_jk2 * value2 + a_jk3 * value3;
3708          }
3709        }
3710#endif
3711      }
3712    }
3713  }
3714}
3715/* Uses factorization to solve. */
3716void ClpCholeskyBase::solve(CoinWorkDouble *region)
3717{
3718  if (!whichDense_) {
3719    solve(region, 3);
3720  } else {
3721    // dense columns
3722    int i;
3723    solve(region, 1);
3724    // do change;
3725    int numberDense = dense_->numberRows();
3726    CoinWorkDouble *change = new CoinWorkDouble[numberDense];
3727    for (i = 0; i < numberDense; i++) {
3728      const longDouble *a = denseColumn_ + i * numberRows_;
3729      CoinWorkDouble value = 0.0;
3730      for (int iRow = 0; iRow < numberRows_; iRow++)
3731        value += a[iRow] * region[iRow];
3732      change[i] = value;
3733    }
3734    // solve
3735    dense_->solve(change);
3736    for (i = 0; i < numberDense; i++) {
3737      const longDouble *a = denseColumn_ + i * numberRows_;
3738      CoinWorkDouble value = change[i];
3739      for (int iRow = 0; iRow < numberRows_; iRow++)
3740        region[iRow] -= value * a[iRow];
3741    }
3742    delete[] change;
3743    // and finish off
3744    solve(region, 2);
3745  }
3746}
3747/* solve - 1 just first half, 2 just second half - 3 both.
3748   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
3749*/
3750void ClpCholeskyBase::solve(CoinWorkDouble *region, int type)
3751{
3752#ifdef CLP_DEBUG
3753  double *regionX = NULL;
3754  if (sizeof(CoinWorkDouble) != sizeof(double) && type == 3) {
3755    regionX = ClpCopyOfArray(region, numberRows_);
3756  }
3757#endif
3758  CoinWorkDouble *work = reinterpret_cast< CoinWorkDouble * >(workDouble_);
3759  int i;
3760  int j;
3761  for (i = 0; i < numberRows_; i++) {
3762    int iRow = permute_[i];
3763    work[i] = region[iRow];
3764  }
3765  switch (type) {
3766  case 1:
3767    for (i = 0; i < numberRows_; i++) {
3768      CoinWorkDouble value = work[i];
3769      int offset = indexStart_[i] - choleskyStart_[i];
3770      for (j = choleskyStart_[i]; j < choleskyStart_[i + 1]; j++) {
3771        int iRow = choleskyRow_[j + offset];
3772        work[iRow] -= sparseFactor_[j] * value;
3773      }
3774    }
3775    for (i = 0; i < numberRows_; i++) {
3776      int iRow = permute_[i];
3777      region[iRow] = work[i] * diagonal_[i];
3778    }
3779    break;
3780  case 2:
3781    for (i = numberRows_ - 1; i >= 0; i--) {
3782      int offset = indexStart_[i] - choleskyStart_[i];
3783      CoinWorkDouble value = work[i] * diagonal_[i];
3784      for (j = choleskyStart_[i]; j < choleskyStart_[i + 1]; j++) {
3785        int iRow = choleskyRow_[j + offset];
3786        value -= sparseFactor_[j] * work[iRow];
3787      }
3788      work[i] = value;
3789      int iRow = permute_[i];
3790      region[iRow] = value;
3791    }
3792    break;
3793  case 3:
3794    for (i = 0; i < firstDense_; i++) {
3795      int offset = indexStart_[i] - choleskyStart_[i];
3796      CoinWorkDouble value = work[i];
3797      for (j = choleskyStart_[i]; j < choleskyStart_[i + 1]; j++) {
3798        int iRow = choleskyRow_[j + offset];
3799        work[iRow] -= sparseFactor_[j] * value;
3800      }
3801    }
3802    if (firstDense_ < numberRows_) {
3803      // do dense
3804      ClpCholeskyDense dense;
3805      // just borrow space
3806      int nDense = numberRows_ - firstDense_;
3807      dense.reserveSpace(this, nDense);
3808      dense.solve(work + firstDense_);
3809      for (i = numberRows_ - 1; i >= firstDense_; i--) {
3810        CoinWorkDouble value = work[i];
3811        int iRow = permute_[i];
3812        region[iRow] = value;
3813      }
3814    }
3815    for (i = firstDense_ - 1; i >= 0; i--) {
3816      int offset = indexStart_[i] - choleskyStart_[i];
3817      CoinWorkDouble value = work[i] * diagonal_[i];
3818      for (j = choleskyStart_[i]; j < choleskyStart_[i + 1]; j++) {
3819        int iRow = choleskyRow_[j + offset];
3820        value -= sparseFactor_[j] * work[iRow];
3821      }
3822      work[i] = value;
3823      int iRow = permute_[i];
3824      region[iRow] = value;
3825    }
3826    break;
3827  }
3828#ifdef CLP_DEBUG
3829  if (regionX) {
3830    longDouble *work = workDouble_;
3831    int i;
3832    int j;
3833    double largestO = 0.0;
3834    for (i = 0; i < numberRows_; i++) {
3835      largestO = CoinMax(largestO, CoinAbs(regionX[i]));
3836    }
3837    for (i = 0; i < numberRows_; i++) {
3838      int iRow = permute_[i];
3839      work[i] = regionX[iRow];
3840    }
3841    for (i = 0; i < firstDense_; i++) {
3842      int offset = indexStart_[i] - choleskyStart_[i];
3843      CoinWorkDouble value = work[i];
3844      for (j = choleskyStart_[i]; j < choleskyStart_[i + 1]; j++) {
3845        int iRow = choleskyRow_[j + offset];
3846        work[iRow] -= sparseFactor_[j] * value;
3847      }
3848    }
3849    if (firstDense_ < numberRows_) {
3850      // do dense
3851      ClpCholeskyDense dense;
3852      // just borrow space
3853      int nDense = numberRows_ - firstDense_;
3854      dense.reserveSpace(this, nDense);
3855      dense.solve(work + firstDense_);
3856      for (i = numberRows_ - 1; i >= firstDense_; i--) {
3857        CoinWorkDouble value = work[i];
3858        int iRow = permute_[i];
3859        regionX[iRow] = value;
3860      }
3861    }
3862    for (i = firstDense_ - 1; i >= 0; i--) {
3863      int offset = indexStart_[i] - choleskyStart_[i];
3864      CoinWorkDouble value = work[i] * diagonal_[i];
3865      for (j = choleskyStart_[i]; j < choleskyStart_[i + 1]; j++) {
3866        int iRow = choleskyRow_[j + offset];
3867        value -= sparseFactor_[j] * work[iRow];
3868      }
3869      work[i] = value;
3870      int iRow = permute_[i];
3871      regionX[iRow] = value;
3872    }
3873    double largest = 0.0;
3874    double largestV = 0.0;
3875    for (i = 0; i < numberRows_; i++) {
3876      largest = CoinMax(largest, CoinAbs(region[i] - regionX[i]));
3877      largestV = CoinMax(largestV, CoinAbs(region[i]));
3878    }
3879    printf("largest difference %g, largest %g, largest original %g\n",
3880      largest, largestV, largestO);
3881    delete[] regionX;
3882  }
3883#endif
3884}
3885#if 0 //CLP_LONG_CHOLESKY
3886/* Uses factorization to solve. */
3887void
3888ClpCholeskyBase::solve (CoinWorkDouble * region)
3889{
3890     assert (!whichDense_) ;
3891     CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
3892     int i;
3893     int j;
3894     for (i = 0; i < numberRows_; i++) {
3895          int iRow = permute_[i];
3896          work[i] = region[iRow];
3897     }
3898     for (i = 0; i < firstDense_; i++) {
3899          int offset = indexStart_[i] - choleskyStart_[i];
3900          CoinWorkDouble value = work[i];
3901          for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
3902               int iRow = choleskyRow_[j+offset];
3903               work[iRow] -= sparseFactor_[j] * value;
3904          }
3905     }
3906     if (firstDense_ < numberRows_) {
3907          // do dense
3908          ClpCholeskyDense dense;
3909          // just borrow space
3910          int nDense = numberRows_ - firstDense_;
3911          dense.reserveSpace(this, nDense);
3912          dense.solve(work + firstDense_);
3913          for (i = numberRows_ - 1; i >= firstDense_; i--) {
3914               CoinWorkDouble value = work[i];
3915               int iRow = permute_[i];
3916               region[iRow] = value;
3917          }
3918     }
3919     for (i = firstDense_ - 1; i >= 0; i--) {
3920          int offset = indexStart_[i] - choleskyStart_[i];
3921          CoinWorkDouble value = work[i] * diagonal_[i];
3922          for (j = choleskyStart_[i]; j < choleskyStart_[i+1]; j++) {
3923               int iRow = choleskyRow_[j+offset];
3924               value -= sparseFactor_[j] * work[iRow];
3925          }
3926          work[i] = value;
3927          int iRow = permute_[i];
3928          region[iRow] = value;
3929     }
3930}
3931#endif
3932
3933/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
3934*/
Note: See TracBrowser for help on using the repository browser.