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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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