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

Last change on this file since 1402 was 1402, checked in by forrest, 10 years ago

get rid of compiler warnings

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