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

Last change on this file since 2030 was 2030, checked in by forrest, 6 years ago

fix some ampl stuff, add ClpSolver? and a few fixes

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