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

Last change on this file since 1368 was 1368, checked in by forrest, 11 years ago

changes for cholesky including code from Anshul Gupta

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