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

Last change on this file since 1878 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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