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

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

try and improve stability - also long double interior point

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 83.0 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <iostream>
7
8#include "ClpCholeskyBase.hpp"
9#include "ClpInterior.hpp"
10#include "ClpHelperFunctions.hpp"
11#include "CoinHelperFunctions.hpp"
12#include "CoinSort.hpp"
13#include "ClpCholeskyDense.hpp"
14#include "ClpMessage.hpp"
15#include "ClpQuadraticObjective.hpp"
16
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24ClpCholeskyBase::ClpCholeskyBase (int denseThreshold) :
25  type_(0),
26  doKKT_(false),
27  goDense_(0.7),
28  choleskyCondition_(0.0),
29  model_(NULL),
30  numberTrials_(),
31  numberRows_(0),
32  status_(0),
33  rowsDropped_(NULL),
34  permuteInverse_(NULL),
35  permute_(NULL),
36  numberRowsDropped_(0),
37  sparseFactor_(NULL),
38  choleskyStart_(NULL),
39  choleskyRow_(NULL),
40  indexStart_(NULL),
41  diagonal_(NULL),
42  workDouble_(NULL),
43  link_(NULL),
44  workInteger_(NULL),
45  clique_(NULL),
46  sizeFactor_(0),
47  sizeIndex_(0),
48  firstDense_(0),
49  rowCopy_(NULL),
50  whichDense_(NULL),
51  denseColumn_(NULL),
52  dense_(NULL),
53  denseThreshold_(denseThreshold)
54{
55  memset(integerParameters_,0,64*sizeof(int));
56  memset(doubleParameters_,0,64*sizeof(double));
57}
58
59//-------------------------------------------------------------------
60// Copy constructor
61//-------------------------------------------------------------------
62ClpCholeskyBase::ClpCholeskyBase (const ClpCholeskyBase & rhs) :
63  type_(rhs.type_),
64  doKKT_(rhs.doKKT_),
65  goDense_(rhs.goDense_),
66  choleskyCondition_(rhs.choleskyCondition_),
67  model_(rhs.model_),
68  numberTrials_(rhs.numberTrials_),
69  numberRows_(rhs.numberRows_),
70  status_(rhs.status_),
71  numberRowsDropped_(rhs.numberRowsDropped_)
72{ 
73  rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_,numberRows_);
74  permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_,numberRows_);
75  permute_ = ClpCopyOfArray(rhs.permute_,numberRows_);
76  sizeFactor_=rhs.sizeFactor_;
77  sizeIndex_ = rhs.sizeIndex_;
78  firstDense_ = rhs.firstDense_;
79  sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
80  choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
81  indexStart_ = ClpCopyOfArray(rhs.indexStart_,numberRows_);
82  choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
83  diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
84#if CLP_LONG_CHOLESKY!=1
85  workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
86#else
87  // actually long double
88  workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
89#endif
90  link_ = ClpCopyOfArray(rhs.link_,numberRows_);
91  workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
92  clique_ = ClpCopyOfArray(rhs.clique_,numberRows_);
93  CoinMemcpyN(rhs.integerParameters_,64,integerParameters_);
94  CoinMemcpyN(rhs.doubleParameters_,64,doubleParameters_);
95  rowCopy_ = rhs.rowCopy_->clone();
96  whichDense_ = NULL;
97  denseColumn_=NULL;
98  dense_=NULL;
99  denseThreshold_ = rhs.denseThreshold_;
100}
101
102//-------------------------------------------------------------------
103// Destructor
104//-------------------------------------------------------------------
105ClpCholeskyBase::~ClpCholeskyBase ()
106{
107  delete [] rowsDropped_;
108  delete [] permuteInverse_;
109  delete [] permute_;
110  delete [] sparseFactor_;
111  delete [] choleskyStart_;
112  delete [] choleskyRow_;
113  delete [] indexStart_;
114  delete [] diagonal_;
115  delete [] workDouble_;
116  delete [] link_;
117  delete [] workInteger_;
118  delete [] clique_;
119  delete rowCopy_;
120  delete [] whichDense_;
121  delete [] denseColumn_;
122  delete dense_;
123}
124
125//----------------------------------------------------------------
126// Assignment operator
127//-------------------------------------------------------------------
128ClpCholeskyBase &
129ClpCholeskyBase::operator=(const ClpCholeskyBase& rhs)
130{
131  if (this != &rhs) {
132    type_ = rhs.type_;
133    doKKT_ = rhs.doKKT_;
134    goDense_ = rhs.goDense_;
135    choleskyCondition_ = rhs.choleskyCondition_;
136    model_ = rhs.model_;
137    numberTrials_ = rhs.numberTrials_;
138    numberRows_ = rhs.numberRows_;
139    status_ = rhs.status_;
140    numberRowsDropped_ = rhs.numberRowsDropped_;
141    delete [] rowsDropped_;
142    delete [] permuteInverse_;
143    delete [] permute_;
144    delete [] sparseFactor_;
145    delete [] choleskyStart_;
146    delete [] choleskyRow_;
147    delete [] indexStart_;
148    delete [] diagonal_;
149    delete [] workDouble_;
150    delete [] link_;
151    delete [] workInteger_;
152    delete [] clique_;
153    delete rowCopy_;
154    delete [] whichDense_;
155    delete [] denseColumn_;
156    delete dense_;
157    rowsDropped_ = ClpCopyOfArray(rhs.rowsDropped_,numberRows_);
158    permuteInverse_ = ClpCopyOfArray(rhs.permuteInverse_,numberRows_);
159    permute_ = ClpCopyOfArray(rhs.permute_,numberRows_);
160    sizeFactor_=rhs.sizeFactor_;
161    sizeIndex_ = rhs.sizeIndex_;
162    firstDense_ = rhs.firstDense_;
163    sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
164    choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
165    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
166    indexStart_ = ClpCopyOfArray(rhs.indexStart_,numberRows_);
167    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,sizeIndex_);
168    diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_);
169#if CLP_LONG_CHOLESKY!=1
170    workDouble_ = ClpCopyOfArray(rhs.workDouble_,numberRows_);
171#else
172    // actually long double
173    workDouble_ = reinterpret_cast<double *> (ClpCopyOfArray(reinterpret_cast<CoinWorkDouble *> (rhs.workDouble_),numberRows_));
174#endif
175    link_ = ClpCopyOfArray(rhs.link_,numberRows_);
176    workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
177    clique_ = ClpCopyOfArray(rhs.clique_,numberRows_);
178    delete rowCopy_;
179    rowCopy_ = rhs.rowCopy_->clone();
180    whichDense_ = NULL;
181    denseColumn_=NULL;
182    dense_=NULL;
183    denseThreshold_ = rhs.denseThreshold_;
184  }
185  return *this;
186}
187// reset numberRowsDropped and rowsDropped.
188void 
189ClpCholeskyBase::resetRowsDropped()
190{
191  numberRowsDropped_=0;
192  memset(rowsDropped_,0,numberRows_);
193}
194/* Uses factorization to solve. - given as if KKT.
195   region1 is rows+columns, region2 is rows */
196void 
197ClpCholeskyBase::solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
198                           CoinWorkDouble diagonalScaleFactor)
199{
200  if (!doKKT_) {
201    int iColumn;
202    int numberColumns = model_->numberColumns();
203    int numberTotal = numberRows_+numberColumns;
204    CoinWorkDouble * region1Save = new CoinWorkDouble[numberTotal];
205    for (iColumn=0;iColumn<numberTotal;iColumn++) {
206      region1[iColumn] *= diagonal[iColumn];
207      region1Save[iColumn]=region1[iColumn];
208    }
209    multiplyAdd(region1+numberColumns,numberRows_,-1.0,region2,1.0);
210    model_->clpMatrix()->times(1.0,region1,region2);
211    CoinWorkDouble maximumRHS = maximumAbsElement(region2,numberRows_);
212    CoinWorkDouble scale=1.0;
213    CoinWorkDouble unscale=1.0;
214    if (maximumRHS>1.0e-30) {
215      if (maximumRHS<=0.5) {
216        CoinWorkDouble factor=2.0;
217        while (maximumRHS<=0.5) {
218          maximumRHS*=factor;
219          scale*=factor;
220        } /* endwhile */
221      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
222        CoinWorkDouble factor=0.5;
223        while (maximumRHS>=2.0) {
224          maximumRHS*=factor;
225          scale*=factor;
226        } /* endwhile */
227      } 
228      unscale=diagonalScaleFactor/scale;
229    } else {
230      //effectively zero
231      scale=0.0;
232      unscale=0.0;
233    } 
234    multiplyAdd(NULL,numberRows_,0.0,region2,scale);
235    solve(region2);
236    multiplyAdd(NULL,numberRows_,0.0,region2,unscale);
237    multiplyAdd(region2,numberRows_,-1.0,region1+numberColumns,0.0);
238    CoinZeroN(region1,numberColumns);
239    model_->clpMatrix()->transposeTimes(1.0,region2,region1);
240    for (iColumn=0;iColumn<numberTotal;iColumn++)
241      region1[iColumn] = region1[iColumn]*diagonal[iColumn]-region1Save[iColumn];
242    delete [] region1Save;
243  } else {
244    // KKT
245    int numberRowsModel = model_->numberRows();
246    int numberColumns = model_->numberColumns();
247    int numberTotal = numberColumns + numberRowsModel;
248    CoinWorkDouble * array = new CoinWorkDouble [numberRows_];
249    CoinMemcpyN(region1,numberTotal,array);
250    CoinMemcpyN(region2,numberRowsModel,array+numberTotal);
251    assert (numberRows_>=numberRowsModel+numberTotal);
252    solve(array);
253    int iRow;
254    for (iRow=0;iRow<numberTotal;iRow++) { 
255      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
256        printf("row region1 %d dropped %g\n",iRow,array[iRow]);
257      }
258    }
259    for (;iRow<numberRows_;iRow++) {
260      if (rowsDropped_[iRow]&&CoinAbs(array[iRow])>1.0e-8) {
261        printf("row region2 %d dropped %g\n",iRow,array[iRow]);
262      }
263    }
264    CoinMemcpyN(array+numberTotal,numberRowsModel,region2);
265    CoinMemcpyN(array,numberTotal,region1);
266    delete [] array;
267  }
268}
269//-------------------------------------------------------------------
270// Clone
271//-------------------------------------------------------------------
272ClpCholeskyBase * ClpCholeskyBase::clone() const
273{
274  return new ClpCholeskyBase(*this);
275}
276// Forms ADAT - returns nonzero if not enough memory
277int 
278ClpCholeskyBase::preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT)
279{
280  delete rowCopy_;
281  rowCopy_ = model_->clpMatrix()->reverseOrderedCopy();
282  if (!doKKT) {
283    numberRows_ = model_->numberRows();
284    rowsDropped_ = new char [numberRows_];
285    memset(rowsDropped_,0,numberRows_);
286    numberRowsDropped_=0;
287    // Space for starts
288    choleskyStart_ = new CoinBigIndex[numberRows_+1];
289    const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
290    const int * columnLength = model_->clpMatrix()->getVectorLengths();
291    const int * row = model_->clpMatrix()->getIndices();
292    const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
293    const int * rowLength = rowCopy_->getVectorLengths();
294    const int * column = rowCopy_->getIndices();
295    // We need two arrays for counts
296    int * which = new int [numberRows_];
297    int * used = new int[numberRows_+1];
298    CoinZeroN(used,numberRows_);
299    int iRow;
300    sizeFactor_=0;
301    int numberColumns = model_->numberColumns();
302    int numberDense=0;
303    //denseThreshold_=3;
304    if (denseThreshold_>0) {
305      delete [] whichDense_;
306      delete [] denseColumn_;
307      delete dense_;
308      whichDense_ = new char[numberColumns];
309      int iColumn;
310      used[numberRows_]=0;
311      for (iColumn=0;iColumn<numberColumns;iColumn++) {
312        int length = columnLength[iColumn];
313        used[length] += 1;
314      }
315      int nLong=0;
316      int stop = CoinMax(denseThreshold_/2,100);
317      for (iRow=numberRows_;iRow>=stop;iRow--) {
318        if (used[iRow]) 
319          printf("%d columns are of length %d\n",used[iRow],iRow);
320        nLong += used[iRow];
321        if (nLong>50||nLong>(numberColumns>>2))
322          break;
323      }
324      CoinZeroN(used,numberRows_);
325      for (iColumn=0;iColumn<numberColumns;iColumn++) {
326        if (columnLength[iColumn]<denseThreshold_) {
327          whichDense_[iColumn]=0;
328        } else {
329          whichDense_[iColumn]=1;
330          numberDense++;
331        }
332      }
333      if (!numberDense||numberDense>100) {
334        // free
335        delete [] whichDense_;
336        whichDense_=NULL;
337        denseColumn_=NULL;
338        dense_=NULL;
339      } else {
340        // space for dense columns
341        denseColumn_ = new longDouble [numberDense*numberRows_];
342        // dense cholesky
343        dense_ = new ClpCholeskyDense();
344        dense_->reserveSpace(NULL,numberDense);
345        printf("Taking %d columns as dense\n",numberDense);
346      }
347    }
348    int offset = includeDiagonal ? 0 : 1;
349    if (lowerTriangular)
350      offset=-offset;
351    for (iRow=0;iRow<numberRows_;iRow++) {
352      int number=0;
353      // make sure diagonal exists if includeDiagonal
354      if (!offset) {
355        which[0]=iRow;
356        used[iRow]=1;
357        number=1;
358      }
359      CoinBigIndex startRow=rowStart[iRow];
360      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
361      if (lowerTriangular) {
362        for (CoinBigIndex k=startRow;k<endRow;k++) {
363          int iColumn=column[k];
364          if (!whichDense_||!whichDense_[iColumn]) {
365            CoinBigIndex start=columnStart[iColumn];
366            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
367            for (CoinBigIndex j=start;j<end;j++) {
368              int jRow=row[j];
369              if (jRow<=iRow+offset) {
370                if (!used[jRow]) {
371                  used[jRow]=1;
372                  which[number++]=jRow;
373                }
374              }
375            }
376          }
377        }
378      } else {
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      }
396      sizeFactor_ += number;
397      int j;
398      for (j=0;j<number;j++)
399        used[which[j]]=0;
400    }
401    delete [] which;
402    // Now we have size - create arrays and fill in
403    try { 
404      choleskyRow_ = new int [sizeFactor_];
405    }
406    catch (...) {
407      // no memory
408      delete [] choleskyStart_;
409      choleskyStart_=NULL;
410      return -1;
411    }
412    sizeFactor_=0;
413    which = choleskyRow_;
414    for (iRow=0;iRow<numberRows_;iRow++) {
415      int number=0;
416      // make sure diagonal exists if includeDiagonal
417      if (!offset) {
418        which[0]=iRow;
419        used[iRow]=1;
420        number=1;
421      }
422      choleskyStart_[iRow]=sizeFactor_;
423      CoinBigIndex startRow=rowStart[iRow];
424      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
425      if (lowerTriangular) {
426        for (CoinBigIndex k=startRow;k<endRow;k++) {
427          int iColumn=column[k];
428          if (!whichDense_||!whichDense_[iColumn]) {
429            CoinBigIndex start=columnStart[iColumn];
430            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
431            for (CoinBigIndex j=start;j<end;j++) {
432              int jRow=row[j];
433              if (jRow<=iRow+offset) {
434                if (!used[jRow]) {
435                  used[jRow]=1;
436                  which[number++]=jRow;
437                }
438              }
439            }
440          }
441        }
442      } else {
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      }
460      sizeFactor_ += number;
461      int j;
462      for (j=0;j<number;j++)
463        used[which[j]]=0;
464      // Sort
465      std::sort(which,which+number);
466      // move which on
467      which += number;
468    }
469    choleskyStart_[numberRows_]=sizeFactor_;
470    delete [] used;
471    return 0;
472  } else {
473    int numberRowsModel = model_->numberRows();
474    int numberColumns = model_->numberColumns();
475    int numberTotal = numberColumns + numberRowsModel;
476    numberRows_ = 2*numberRowsModel+numberColumns;
477    rowsDropped_ = new char [numberRows_];
478    memset(rowsDropped_,0,numberRows_);
479    numberRowsDropped_=0;
480    CoinPackedMatrix * quadratic = NULL;
481    ClpQuadraticObjective * quadraticObj = 
482      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
483    if (quadraticObj) 
484      quadratic = quadraticObj->quadraticObjective();
485    int numberElements = model_->clpMatrix()->getNumElements();
486    numberElements = numberElements+2*numberRowsModel+numberTotal;
487    if (quadratic)
488      numberElements += quadratic->getNumElements(); 
489    // Space for starts
490    choleskyStart_ = new CoinBigIndex[numberRows_+1];
491    const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
492    const int * columnLength = model_->clpMatrix()->getVectorLengths();
493    const int * row = model_->clpMatrix()->getIndices();
494    //const double * element = model_->clpMatrix()->getElements();
495    // Now we have size - create arrays and fill in
496    try { 
497      choleskyRow_ = new int [numberElements];
498    }
499    catch (...) {
500      // no memory
501      delete [] choleskyStart_;
502      choleskyStart_=NULL;
503      return -1;
504    }
505    int iRow,iColumn;
506 
507    sizeFactor_=0;
508    // matrix
509    if (lowerTriangular) {
510      if (!quadratic) {
511        for (iColumn=0;iColumn<numberColumns;iColumn++) {
512          choleskyStart_[iColumn]=sizeFactor_;
513          choleskyRow_[sizeFactor_++]=iColumn;
514          CoinBigIndex start=columnStart[iColumn];
515          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
516          if (!includeDiagonal)
517            start++;
518          for (CoinBigIndex j=start;j<end;j++) {
519            choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
520          }
521        }
522      } else {
523        // Quadratic
524        const int * columnQuadratic = quadratic->getIndices();
525        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
526        const int * columnQuadraticLength = quadratic->getVectorLengths();
527        //const double * quadraticElement = quadratic->getElements();
528        for (iColumn=0;iColumn<numberColumns;iColumn++) {
529          choleskyStart_[iColumn]=sizeFactor_;
530          if (includeDiagonal) 
531            choleskyRow_[sizeFactor_++]=iColumn;
532    CoinBigIndex j;
533          for ( j=columnQuadraticStart[iColumn];
534               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
535            int jColumn = columnQuadratic[j];
536            if (jColumn>iColumn)
537              choleskyRow_[sizeFactor_++]=jColumn;
538          }
539          CoinBigIndex start=columnStart[iColumn];
540          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
541          for ( j=start;j<end;j++) {
542            choleskyRow_[sizeFactor_++]=row[j]+numberTotal;
543          }
544        }
545      }
546      // slacks
547      for (;iColumn<numberTotal;iColumn++) {
548        choleskyStart_[iColumn]=sizeFactor_;
549        if (includeDiagonal) 
550          choleskyRow_[sizeFactor_++]=iColumn;
551        choleskyRow_[sizeFactor_++]=iColumn-numberColumns+numberTotal;
552      }
553      // Transpose - nonzero diagonal (may regularize)
554      for (iRow=0;iRow<numberRowsModel;iRow++) {
555        choleskyStart_[iRow+numberTotal]=sizeFactor_;
556        // diagonal
557        if (includeDiagonal) 
558          choleskyRow_[sizeFactor_++]=iRow+numberTotal;
559      }
560      choleskyStart_[numberRows_]=sizeFactor_;
561    } else {
562      // transpose
563      ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
564      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
565      const int * rowLength = rowCopy->getVectorLengths();
566      const int * column = rowCopy->getIndices();
567      if (!quadratic) {
568        for (iColumn=0;iColumn<numberColumns;iColumn++) {
569          choleskyStart_[iColumn]=sizeFactor_;
570          if (includeDiagonal) 
571            choleskyRow_[sizeFactor_++]=iColumn;
572        }
573      } else {
574        // Quadratic
575        // transpose
576        CoinPackedMatrix quadraticT;
577        quadraticT.reverseOrderedCopyOf(*quadratic);
578        const int * columnQuadratic = quadraticT.getIndices();
579        const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
580        const int * columnQuadraticLength = quadraticT.getVectorLengths();
581        //const double * quadraticElement = quadraticT.getElements();
582        for (iColumn=0;iColumn<numberColumns;iColumn++) {
583          choleskyStart_[iColumn]=sizeFactor_;
584          for (CoinBigIndex j=columnQuadraticStart[iColumn];
585               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
586            int jColumn = columnQuadratic[j];
587            if (jColumn<iColumn)
588              choleskyRow_[sizeFactor_++]=jColumn;
589          }
590          if (includeDiagonal) 
591            choleskyRow_[sizeFactor_++]=iColumn;
592        }
593      }
594      int iRow;
595      // slacks
596      for (iRow=0;iRow<numberRowsModel;iRow++) {
597        choleskyStart_[iRow+numberColumns]=sizeFactor_;
598        if (includeDiagonal) 
599          choleskyRow_[sizeFactor_++]=iRow+numberColumns;
600      }
601      for (iRow=0;iRow<numberRowsModel;iRow++) {
602        choleskyStart_[iRow+numberTotal]=sizeFactor_;
603        CoinBigIndex start=rowStart[iRow];
604        CoinBigIndex end=rowStart[iRow]+rowLength[iRow];
605        for (CoinBigIndex j=start;j<end;j++) {
606          choleskyRow_[sizeFactor_++]=column[j];
607        }
608        // slack
609        choleskyRow_[sizeFactor_++]=numberColumns+iRow;
610        if (includeDiagonal)
611          choleskyRow_[sizeFactor_++]=iRow+numberTotal;
612      }
613      choleskyStart_[numberRows_]=sizeFactor_;
614    }
615  }
616  return 0;
617}
618/* Orders rows and saves pointer to matrix.and model */
619int 
620ClpCholeskyBase::order(ClpInterior * model) 
621{
622  model_=model;
623  int numberRowsModel = model_->numberRows();
624  int numberColumns = model_->numberColumns();
625  int numberTotal = numberColumns + numberRowsModel;
626  CoinPackedMatrix * quadratic = NULL;
627  ClpQuadraticObjective * quadraticObj = 
628    (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
629  if (quadraticObj) 
630    quadratic = quadraticObj->quadraticObjective();
631  if (!doKKT_) {
632    numberRows_ = model->numberRows();
633  } else {
634    numberRows_ = 2*numberRowsModel+numberColumns;
635  }
636  rowsDropped_ = new char [numberRows_];
637  numberRowsDropped_=0;
638  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
639  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
640  const int * columnLength = model_->clpMatrix()->getVectorLengths();
641  const int * row = model_->clpMatrix()->getIndices();
642  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
643  const int * rowLength = rowCopy_->getVectorLengths();
644  const int * column = rowCopy_->getIndices();
645  // We need arrays for counts
646  int * which = new int [numberRows_];
647  int * used = new int[numberRows_+1];
648  int * count = new int[numberRows_];
649  CoinZeroN(count,numberRows_);
650  CoinZeroN(used,numberRows_);
651  int iRow;
652  sizeFactor_=0;
653  permute_ = new int[numberRows_];
654  for (iRow=0;iRow<numberRows_;iRow++) 
655    permute_[iRow]=iRow; 
656  if (!doKKT_) {
657    int numberDense=0;
658    if (denseThreshold_>0) {
659      delete [] whichDense_;
660      delete [] denseColumn_;
661      delete dense_;
662      whichDense_ = new char[numberColumns];
663      int iColumn;
664      used[numberRows_]=0;
665      for (iColumn=0;iColumn<numberColumns;iColumn++) {
666        int length = columnLength[iColumn];
667        used[length] += 1;
668      }
669      int nLong=0;
670      int stop = CoinMax(denseThreshold_/2,100);
671      for (iRow=numberRows_;iRow>=stop;iRow--) {
672        if (used[iRow]) 
673          printf("%d columns are of length %d\n",used[iRow],iRow);
674        nLong += used[iRow];
675        if (nLong>50||nLong>(numberColumns>>2))
676          break;
677      }
678      CoinZeroN(used,numberRows_);
679      for (iColumn=0;iColumn<numberColumns;iColumn++) {
680        if (columnLength[iColumn]<denseThreshold_) {
681          whichDense_[iColumn]=0;
682        } else {
683          whichDense_[iColumn]=1;
684          numberDense++;
685        }
686      }
687      if (!numberDense||numberDense>100) {
688        // free
689        delete [] whichDense_;
690        whichDense_=NULL;
691        denseColumn_=NULL;
692        dense_=NULL;
693      } else {
694        // space for dense columns
695        denseColumn_ = new longDouble [numberDense*numberRows_];
696        // dense cholesky
697        dense_ = new ClpCholeskyDense();
698        dense_->reserveSpace(NULL,numberDense);
699        printf("Taking %d columns as dense\n",numberDense);
700      }
701    }
702    /*
703       Get row counts and size
704    */
705    for (iRow=0;iRow<numberRows_;iRow++) {
706      int number=1;
707      // make sure diagonal exists
708      which[0]=iRow;
709      used[iRow]=1;
710      CoinBigIndex startRow=rowStart[iRow];
711      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
712      for (CoinBigIndex k=startRow;k<endRow;k++) {
713        int iColumn=column[k];
714        if (!whichDense_||!whichDense_[iColumn]) {
715          CoinBigIndex start=columnStart[iColumn];
716          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
717          for (CoinBigIndex j=start;j<end;j++) {
718            int jRow=row[j];
719            if (jRow<iRow) {
720              if (!used[jRow]) {
721                used[jRow]=1;
722                which[number++]=jRow;
723                count[jRow]++;
724              }
725            }
726          }
727        }
728      }
729      sizeFactor_ += number;
730      count[iRow]+=number;
731      int j;
732      for (j=0;j<number;j++)
733        used[which[j]]=0;
734    }
735    CoinSort_2(count,count+numberRows_,permute_);
736  } else {
737    // KKT
738    int numberElements = model_->clpMatrix()->getNumElements();
739    numberElements = numberElements+2*numberRowsModel+numberTotal;
740    if (quadratic)
741      numberElements += quadratic->getNumElements(); 
742    // off diagonal
743    numberElements -= numberRows_;
744    sizeFactor_=numberElements;
745    // If we sort we need to redo symbolic etc
746  }
747  delete [] which;
748  delete [] used;
749  delete [] count;
750  permuteInverse_ = new int [numberRows_];
751  memset(rowsDropped_,0,numberRows_);
752  for (iRow=0;iRow<numberRows_;iRow++) {
753    //permute_[iRow]=iRow; // force no permute
754    //permute_[iRow]=numberRows_-1-iRow; // force odd permute
755    //permute_[iRow]=(iRow+1)%numberRows_; // force odd permute
756    permuteInverse_[permute_[iRow]]=iRow;
757  }
758  return 0;
759}
760/* Does Symbolic factorization given permutation.
761   This is called immediately after order.  If user provides this then
762   user must provide factorize and solve.  Otherwise the default factorization is used
763   returns non-zero if not enough memory */
764int 
765ClpCholeskyBase::symbolic()
766{
767  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
768  const int * columnLength = model_->clpMatrix()->getVectorLengths();
769  const int * row = model_->clpMatrix()->getIndices();
770  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
771  const int * rowLength = rowCopy_->getVectorLengths();
772  const int * column = rowCopy_->getIndices();
773  int numberRowsModel = model_->numberRows();
774  int numberColumns = model_->numberColumns();
775  int numberTotal = numberColumns + numberRowsModel;
776  CoinPackedMatrix * quadratic = NULL;
777  ClpQuadraticObjective * quadraticObj = 
778    (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
779  if (quadraticObj) 
780    quadratic = quadraticObj->quadraticObjective();
781  // We need an array for counts
782  int * used = new int[numberRows_+1];
783  // If KKT then re-order so negative first
784  if (doKKT_) {
785    int nn=0;
786    int np=0;
787    int iRow;
788    for (iRow=0;iRow<numberRows_;iRow++) {
789      int originalRow = permute_[iRow];
790      if (originalRow<numberTotal)
791        permute_[nn++]=originalRow;
792      else
793        used[np++]=originalRow;
794    }
795    CoinMemcpyN(used,np,permute_+nn);
796    for (iRow=0;iRow<numberRows_;iRow++) 
797      permuteInverse_[permute_[iRow]]=iRow;
798  }
799  CoinZeroN(used,numberRows_);
800  int iRow;
801  int iColumn;
802  bool noMemory=false;
803  CoinBigIndex * Astart = new CoinBigIndex[numberRows_+1];
804  int * Arow=NULL;
805  try { 
806    Arow = new int [sizeFactor_];
807  }
808  catch (...) {
809    // no memory
810    delete [] Astart;
811    return -1;
812  }
813  choleskyStart_ = new int[numberRows_+1];
814  link_ = new int[numberRows_];
815  workInteger_ = new CoinBigIndex[numberRows_];
816  indexStart_ = new CoinBigIndex[numberRows_];
817  clique_ = new int[numberRows_];
818  // Redo so permuted upper triangular
819  sizeFactor_=0;
820  int * which = Arow;
821  if (!doKKT_) {
822    for (iRow=0;iRow<numberRows_;iRow++) {
823      int number=0;
824      int iOriginalRow = permute_[iRow];
825      Astart[iRow]=sizeFactor_;
826      CoinBigIndex startRow=rowStart[iOriginalRow];
827      CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
828      for (CoinBigIndex k=startRow;k<endRow;k++) {
829        int iColumn=column[k];
830        if (!whichDense_||!whichDense_[iColumn]) {
831          CoinBigIndex start=columnStart[iColumn];
832          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
833          for (CoinBigIndex j=start;j<end;j++) {
834            int jRow=row[j];
835            int jNewRow = permuteInverse_[jRow];
836            if (jNewRow<iRow) {
837              if (!used[jNewRow]) {
838                used[jNewRow]=1;
839                which[number++]=jNewRow;
840              }
841            }
842          }
843        }
844      }
845      sizeFactor_ += number;
846      int j;
847      for (j=0;j<number;j++)
848        used[which[j]]=0;
849      // Sort
850      std::sort(which,which+number);
851      // move which on
852      which += number;
853    }
854  } else {
855    // KKT
856    // transpose
857    ClpMatrixBase * rowCopy = model_->clpMatrix()->reverseOrderedCopy();
858    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
859    const int * rowLength = rowCopy->getVectorLengths();
860    const int * column = rowCopy->getIndices();
861    // temp
862    bool permuted=false;
863    for (iRow=0;iRow<numberRows_;iRow++) {
864      if (permute_[iRow]!=iRow) {
865        permuted=true;
866        break;
867      }
868    }
869    if (permuted) {
870      // Need to permute - ugly
871      if (!quadratic) {
872        for (iRow=0;iRow<numberRows_;iRow++) {
873          Astart[iRow]=sizeFactor_;
874          int iOriginalRow = permute_[iRow];
875          if (iOriginalRow<numberColumns) {
876            // A may be upper triangular by mistake
877            iColumn=iOriginalRow;
878            CoinBigIndex start=columnStart[iColumn];
879            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
880            for (CoinBigIndex j=start;j<end;j++) {
881              int kRow = row[j]+numberTotal;
882              kRow=permuteInverse_[kRow];
883              if (kRow<iRow)
884                Arow[sizeFactor_++]=kRow;
885            }
886          } else if (iOriginalRow<numberTotal) {
887            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
888            if (kRow<iRow)
889              Arow[sizeFactor_++]=kRow;
890          } else {
891            int kRow = iOriginalRow-numberTotal;
892            CoinBigIndex start=rowStart[kRow];
893            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
894            for (CoinBigIndex j=start;j<end;j++) {
895              int jRow = column[j];
896              int jNewRow = permuteInverse_[jRow];
897              if (jNewRow<iRow)
898                Arow[sizeFactor_++]=jNewRow;
899            }
900            // slack - should it be permute
901            kRow = permuteInverse_[kRow+numberColumns];
902            if (kRow<iRow)
903              Arow[sizeFactor_++]=kRow;
904          }
905          // Sort
906          std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
907        }
908      } else {
909        // quadratic
910        // transpose
911        CoinPackedMatrix quadraticT;
912        quadraticT.reverseOrderedCopyOf(*quadratic);
913        const int * columnQuadratic = quadraticT.getIndices();
914        const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
915        const int * columnQuadraticLength = quadraticT.getVectorLengths();
916        for (iRow=0;iRow<numberRows_;iRow++) {
917          Astart[iRow]=sizeFactor_;
918          int iOriginalRow = permute_[iRow];
919          if (iOriginalRow<numberColumns) {
920            // Quadratic bit
921            CoinBigIndex j;
922            for ( j=columnQuadraticStart[iOriginalRow];
923                  j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
924              int jColumn = columnQuadratic[j];
925              int jNewColumn = permuteInverse_[jColumn];
926              if (jNewColumn<iRow)
927                Arow[sizeFactor_++]=jNewColumn;
928            }
929            // A may be upper triangular by mistake
930            iColumn=iOriginalRow;
931            CoinBigIndex start=columnStart[iColumn];
932            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
933            for (j=start;j<end;j++) {
934              int kRow = row[j]+numberTotal;
935              kRow=permuteInverse_[kRow];
936              if (kRow<iRow)
937                Arow[sizeFactor_++]=kRow;
938            }
939          } else if (iOriginalRow<numberTotal) {
940            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
941            if (kRow<iRow)
942              Arow[sizeFactor_++]=kRow;
943          } else {
944            int kRow = iOriginalRow-numberTotal;
945            CoinBigIndex start=rowStart[kRow];
946            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
947            for (CoinBigIndex j=start;j<end;j++) {
948              int jRow = column[j];
949              int jNewRow = permuteInverse_[jRow];
950              if (jNewRow<iRow)
951                Arow[sizeFactor_++]=jNewRow;
952            }
953            // slack - should it be permute
954            kRow = permuteInverse_[kRow+numberColumns];
955            if (kRow<iRow)
956              Arow[sizeFactor_++]=kRow;
957          }
958          // Sort
959          std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
960        }
961      }
962    } else {
963      if (!quadratic) {
964        for (iRow=0;iRow<numberRows_;iRow++) {
965          Astart[iRow]=sizeFactor_;
966        }
967      } else {
968        // Quadratic
969        // transpose
970        CoinPackedMatrix quadraticT;
971        quadraticT.reverseOrderedCopyOf(*quadratic);
972        const int * columnQuadratic = quadraticT.getIndices();
973        const CoinBigIndex * columnQuadraticStart = quadraticT.getVectorStarts();
974        const int * columnQuadraticLength = quadraticT.getVectorLengths();
975        //const double * quadraticElement = quadraticT.getElements();
976        for (iRow=0;iRow<numberColumns;iRow++) {
977          int iOriginalRow = permute_[iRow];
978          Astart[iRow]=sizeFactor_;
979          for (CoinBigIndex j=columnQuadraticStart[iOriginalRow];
980               j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
981            int jColumn = columnQuadratic[j];
982            int jNewColumn = permuteInverse_[jColumn];
983            if (jNewColumn<iRow)
984              Arow[sizeFactor_++]=jNewColumn;
985          }
986        }
987      }
988      int iRow;
989      // slacks
990      for (iRow=0;iRow<numberRowsModel;iRow++) {
991        Astart[iRow+numberColumns]=sizeFactor_;
992      }
993      for (iRow=0;iRow<numberRowsModel;iRow++) {
994        Astart[iRow+numberTotal]=sizeFactor_;
995        CoinBigIndex start=rowStart[iRow];
996        CoinBigIndex end=rowStart[iRow]+rowLength[iRow];
997        for (CoinBigIndex j=start;j<end;j++) {
998          Arow[sizeFactor_++]=column[j];
999        }
1000        // slack
1001        Arow[sizeFactor_++]=numberColumns+iRow;
1002      }
1003    }
1004    delete rowCopy;
1005  }
1006  Astart[numberRows_]=sizeFactor_;
1007  firstDense_=numberRows_;
1008  symbolic1(Astart,Arow);
1009  // Now fill in indices
1010  try { 
1011    // too big
1012    choleskyRow_ = new int[sizeFactor_];
1013  }
1014  catch (...) {
1015    // no memory
1016    noMemory=true;
1017  } 
1018  double sizeFactor=sizeFactor_;
1019  if (!noMemory) {
1020    // Do lower triangular
1021    sizeFactor_=0;
1022    int * which = Arow;
1023    if (!doKKT_) {
1024      for (iRow=0;iRow<numberRows_;iRow++) {
1025        int number=0;
1026        int iOriginalRow = permute_[iRow];
1027        Astart[iRow]=sizeFactor_;
1028        if (!rowsDropped_[iOriginalRow]) {
1029          CoinBigIndex startRow=rowStart[iOriginalRow];
1030          CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
1031          for (CoinBigIndex k=startRow;k<endRow;k++) {
1032            int iColumn=column[k];
1033            if (!whichDense_||!whichDense_[iColumn]) {
1034              CoinBigIndex start=columnStart[iColumn];
1035              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1036              for (CoinBigIndex j=start;j<end;j++) {
1037                int jRow=row[j];
1038                int jNewRow = permuteInverse_[jRow];
1039                if (jNewRow>iRow&&!rowsDropped_[jRow]) {
1040                  if (!used[jNewRow]) {
1041                    used[jNewRow]=1;
1042                    which[number++]=jNewRow;
1043                  }
1044                }
1045              }
1046            }
1047          }
1048          sizeFactor_ += number;
1049          int j;
1050          for (j=0;j<number;j++)
1051            used[which[j]]=0;
1052          // Sort
1053          std::sort(which,which+number);
1054          // move which on
1055          which += number;
1056        }
1057      }
1058    } else {
1059      // KKT
1060      // temp
1061      bool permuted=false;
1062      for (iRow=0;iRow<numberRows_;iRow++) {
1063        if (permute_[iRow]!=iRow) {
1064          permuted=true;
1065          break;
1066        }
1067      }
1068      // but fake it
1069      for (iRow=0;iRow<numberRows_;iRow++) {
1070        //permute_[iRow]=iRow; // force no permute
1071        //permuteInverse_[permute_[iRow]]=iRow;
1072      }
1073      if (permuted) {
1074        // Need to permute - ugly
1075        if (!quadratic) {
1076          for (iRow=0;iRow<numberRows_;iRow++) {
1077            Astart[iRow]=sizeFactor_;
1078            int iOriginalRow = permute_[iRow];
1079            if (iOriginalRow<numberColumns) {
1080              iColumn=iOriginalRow;
1081              CoinBigIndex start=columnStart[iColumn];
1082              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1083              for (CoinBigIndex j=start;j<end;j++) {
1084                int kRow = row[j]+numberTotal;
1085                kRow=permuteInverse_[kRow];
1086                if (kRow>iRow)
1087                  Arow[sizeFactor_++]=kRow;
1088              }
1089            } else if (iOriginalRow<numberTotal) {
1090              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1091              if (kRow>iRow)
1092                Arow[sizeFactor_++]=kRow;
1093            } else {
1094              int kRow = iOriginalRow-numberTotal;
1095              CoinBigIndex start=rowStart[kRow];
1096              CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1097              for (CoinBigIndex j=start;j<end;j++) {
1098                int jRow = column[j];
1099                int jNewRow = permuteInverse_[jRow];
1100                if (jNewRow>iRow)
1101                  Arow[sizeFactor_++]=jNewRow;
1102              }
1103              // slack - should it be permute
1104              kRow = permuteInverse_[kRow+numberColumns];
1105              if (kRow>iRow)
1106                Arow[sizeFactor_++]=kRow;
1107            }
1108            // Sort
1109            std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
1110          }
1111        } else {
1112          // quadratic
1113          const int * columnQuadratic = quadratic->getIndices();
1114          const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1115          const int * columnQuadraticLength = quadratic->getVectorLengths();
1116          for (iRow=0;iRow<numberRows_;iRow++) {
1117            Astart[iRow]=sizeFactor_;
1118            int iOriginalRow = permute_[iRow];
1119            if (iOriginalRow<numberColumns) {
1120              // Quadratic bit
1121              CoinBigIndex j;
1122              for ( j=columnQuadraticStart[iOriginalRow];
1123                    j<columnQuadraticStart[iOriginalRow]+columnQuadraticLength[iOriginalRow];j++) {
1124                int jColumn = columnQuadratic[j];
1125                int jNewColumn = permuteInverse_[jColumn];
1126                if (jNewColumn>iRow)
1127                  Arow[sizeFactor_++]=jNewColumn;
1128              }
1129              iColumn=iOriginalRow;
1130              CoinBigIndex start=columnStart[iColumn];
1131              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1132              for (j=start;j<end;j++) {
1133                int kRow = row[j]+numberTotal;
1134                kRow=permuteInverse_[kRow];
1135                if (kRow>iRow)
1136                  Arow[sizeFactor_++]=kRow;
1137              }
1138            } else if (iOriginalRow<numberTotal) {
1139              int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1140              if (kRow>iRow)
1141                Arow[sizeFactor_++]=kRow;
1142            } else {
1143              int kRow = iOriginalRow-numberTotal;
1144              CoinBigIndex start=rowStart[kRow];
1145              CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1146              for (CoinBigIndex j=start;j<end;j++) {
1147                int jRow = column[j];
1148                int jNewRow = permuteInverse_[jRow];
1149                if (jNewRow>iRow)
1150                  Arow[sizeFactor_++]=jNewRow;
1151              }
1152              // slack - should it be permute
1153              kRow = permuteInverse_[kRow+numberColumns];
1154              if (kRow>iRow)
1155                Arow[sizeFactor_++]=kRow;
1156            }
1157            // Sort
1158            std::sort(Arow+Astart[iRow],Arow+sizeFactor_);
1159          }
1160        }
1161      } else {
1162        if (!quadratic) {
1163          for (iColumn=0;iColumn<numberColumns;iColumn++) {
1164            Astart[iColumn]=sizeFactor_;
1165            CoinBigIndex start=columnStart[iColumn];
1166            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1167            for (CoinBigIndex j=start;j<end;j++) {
1168              Arow[sizeFactor_++]=row[j]+numberTotal;
1169            }
1170          }
1171        } else {
1172          // Quadratic
1173          const int * columnQuadratic = quadratic->getIndices();
1174          const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1175          const int * columnQuadraticLength = quadratic->getVectorLengths();
1176          //const double * quadraticElement = quadratic->getElements();
1177          for (iColumn=0;iColumn<numberColumns;iColumn++) {
1178            Astart[iColumn]=sizeFactor_;
1179      CoinBigIndex j;
1180            for ( j=columnQuadraticStart[iColumn];
1181                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1182              int jColumn = columnQuadratic[j];
1183              if (jColumn>iColumn)
1184                Arow[sizeFactor_++]=jColumn;
1185            }
1186            CoinBigIndex start=columnStart[iColumn];
1187            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1188            for ( j=start;j<end;j++) {
1189              Arow[sizeFactor_++]=row[j]+numberTotal;
1190            }
1191          }
1192        }
1193        // slacks
1194        for (iRow=0;iRow<numberRowsModel;iRow++) {
1195          Astart[iRow+numberColumns]=sizeFactor_;
1196          Arow[sizeFactor_++]=iRow+numberTotal;
1197        }
1198        // Transpose - nonzero diagonal (may regularize)
1199        for (iRow=0;iRow<numberRowsModel;iRow++) {
1200          Astart[iRow+numberTotal]=sizeFactor_;
1201        }
1202      }
1203      Astart[numberRows_]=sizeFactor_;
1204    }
1205    symbolic2(Astart,Arow);
1206    if (sizeIndex_<sizeFactor_) {
1207      int * indices=NULL;
1208      try { 
1209        indices = new int[sizeIndex_];
1210      }
1211      catch (...) {
1212        // no memory
1213        noMemory=true;
1214      } 
1215      if (!noMemory)  {
1216 CoinMemcpyN(choleskyRow_,sizeIndex_,indices);
1217        delete [] choleskyRow_;
1218        choleskyRow_=indices;
1219      }
1220    }
1221  }
1222  delete [] used;
1223  // Use cholesky regions
1224  delete [] Astart;
1225  delete [] Arow;
1226  double flops=0.0;
1227  for (iRow=0;iRow<numberRows_;iRow++) {
1228    int length = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1229    flops += static_cast<double> (length) * (length + 2.0);
1230  }
1231  if (model_->messageHandler()->logLevel()>0) 
1232    std::cout<<sizeFactor<<" elements in sparse Cholesky, flop count "<<flops<<std::endl;
1233  try { 
1234    sparseFactor_ = new longDouble [sizeFactor_];
1235#if CLP_LONG_CHOLESKY!=1
1236    workDouble_ = new longDouble[numberRows_];
1237#else
1238    // actually long double
1239    workDouble_ = reinterpret_cast<double *> (new CoinWorkDouble[numberRows_]);
1240#endif
1241    diagonal_ = new longDouble[numberRows_];
1242  }
1243  catch (...) {
1244    // no memory
1245    noMemory=true;
1246  }
1247  if (noMemory) {
1248    delete [] choleskyRow_;
1249    choleskyRow_=NULL;
1250    delete [] choleskyStart_;
1251    choleskyStart_=NULL;
1252    delete [] permuteInverse_;
1253    permuteInverse_ = NULL;
1254    delete [] permute_;
1255    permute_ = NULL;
1256    delete [] choleskyStart_;
1257    choleskyStart_ = NULL;
1258    delete [] indexStart_;
1259    indexStart_ = NULL;
1260    delete [] link_;
1261    link_ = NULL;
1262    delete [] workInteger_;
1263    workInteger_ = NULL;
1264    delete [] sparseFactor_;
1265    sparseFactor_ = NULL;
1266    delete [] workDouble_;
1267    workDouble_ = NULL;
1268    delete [] diagonal_;
1269    diagonal_ = NULL;
1270    delete [] clique_;
1271    clique_ = NULL;
1272    return -1;
1273  }
1274  return 0;
1275}
1276int
1277ClpCholeskyBase::symbolic1(const CoinBigIndex * Astart, const int * Arow)
1278{
1279  int * marked = reinterpret_cast<int *> (workInteger_);
1280  int iRow;
1281  // may not need to do this here but makes debugging easier
1282  for (iRow=0;iRow<numberRows_;iRow++) {
1283    marked[iRow]=-1;
1284    link_[iRow]=-1;
1285    choleskyStart_[iRow]=0; // counts
1286  }
1287  for (iRow=0;iRow<numberRows_;iRow++) {
1288    marked[iRow]=iRow;
1289    for (CoinBigIndex j=Astart[iRow];j<Astart[iRow+1];j++) {
1290      int kRow=Arow[j];
1291      while (marked[kRow] != iRow ) {
1292        if (link_[kRow] <0 )
1293          link_[kRow]=iRow;
1294        choleskyStart_[kRow]++;
1295        marked[kRow]=iRow;
1296        kRow = link_[kRow];
1297      }
1298    }
1299  }
1300  sizeFactor_=0;
1301  for (iRow=0;iRow<numberRows_;iRow++) {
1302    int number = choleskyStart_[iRow];
1303    choleskyStart_[iRow]=sizeFactor_;
1304    sizeFactor_ += number;
1305  }
1306  choleskyStart_[numberRows_]=sizeFactor_;
1307  return sizeFactor_;;
1308}
1309void
1310ClpCholeskyBase::symbolic2(const CoinBigIndex * Astart, const int * Arow)
1311{
1312  int * mergeLink = clique_;
1313  int * marker = reinterpret_cast<int *> (workInteger_);
1314  int iRow;
1315  for (iRow=0;iRow<numberRows_;iRow++) {
1316    marker[iRow]=-1;
1317    mergeLink[iRow]=-1;
1318    link_[iRow]=-1; // not needed but makes debugging easier
1319  }
1320  int start=0;
1321  int end=0;
1322  choleskyStart_[0]=0;
1323   
1324  for (iRow=0;iRow<numberRows_;iRow++) {
1325    int nz=0;
1326    int merge = mergeLink[iRow];
1327    bool marked=false;
1328    if (merge<0)
1329      marker[iRow]=iRow;
1330    else
1331      marker[iRow]=merge;
1332    start = end;
1333    int startSub=start;
1334    link_[iRow]=numberRows_;
1335    CoinBigIndex j;
1336    for ( j=Astart[iRow];j<Astart[iRow+1];j++) {
1337      int kRow=Arow[j];
1338      int k=iRow;
1339      int linked = link_[iRow];
1340      while (linked<=kRow) {
1341        k=linked;
1342        linked = link_[k];
1343      }
1344      nz++;
1345      link_[k]=kRow;
1346      link_[kRow]=linked;
1347      if (marker[kRow] != marker[iRow])
1348        marked=true;
1349    }
1350    bool reuse=false;
1351    // Check if we can re-use indices
1352    if (!marked && merge>=0&&mergeLink[merge]<0) {
1353      // can re-use all
1354      startSub = indexStart_[merge]+1;
1355      nz=choleskyStart_[merge+1]-(choleskyStart_[merge]+1);
1356      reuse=true;
1357    } else {
1358      // See if we can re-use any
1359      int k=mergeLink[iRow];
1360      int maxLength=0;
1361      while (k>=0) {
1362        int length = choleskyStart_[k+1]-(choleskyStart_[k]+1);
1363        int start = indexStart_[k]+1;
1364        int stop = start+length;
1365        if (length>maxLength) {
1366          maxLength = length;
1367          startSub = start;
1368        }
1369        int linked = iRow;
1370        for (CoinBigIndex j=start;j<stop;j++) {
1371          int kRow=choleskyRow_[j];
1372          int kk=linked;
1373          linked = link_[kk];
1374          while (linked<kRow) {
1375            kk=linked;
1376            linked = link_[kk];
1377          }
1378          if (linked!=kRow) {
1379            nz++;
1380            link_[kk]=kRow;
1381            link_[kRow]=linked;
1382            linked=kRow;
1383          }
1384        }
1385        k=mergeLink[k];
1386      }
1387      if (nz== maxLength) 
1388        reuse=true; // can re-use
1389    }
1390    //reuse=false; //temp
1391    if (!reuse) {
1392      end += nz;
1393      startSub=start;
1394      int kRow = iRow;
1395      for (int j=start;j<end;j++) {
1396        kRow=link_[kRow];
1397        choleskyRow_[j]=kRow;
1398        assert (kRow<numberRows_);
1399        marker[kRow]=iRow;
1400      }
1401      marker[iRow]=iRow;
1402    }
1403    indexStart_[iRow]=startSub;
1404    choleskyStart_[iRow+1]=choleskyStart_[iRow]+nz;
1405    if (nz>1) {
1406      int kRow = choleskyRow_[startSub];
1407      mergeLink[iRow]=mergeLink[kRow];
1408      mergeLink[kRow]=iRow;
1409    }
1410    // should not be needed
1411    //std::sort(choleskyRow_+indexStart_[iRow]
1412    //      ,choleskyRow_+indexStart_[iRow]+nz);
1413    //#define CLP_DEBUG
1414#ifdef CLP_DEBUG
1415    int last=-1;
1416    for ( j=indexStart_[iRow];j<indexStart_[iRow]+nz;j++) {
1417      int kRow=choleskyRow_[j];
1418      assert (kRow>last);
1419      last=kRow;
1420    }
1421#endif   
1422  }
1423  sizeFactor_ = choleskyStart_[numberRows_];
1424  sizeIndex_ = start;
1425  // find dense segment here
1426  int numberleft=numberRows_;
1427  for (iRow=0;iRow<numberRows_;iRow++) {
1428    CoinBigIndex left=sizeFactor_-choleskyStart_[iRow];
1429    double n=numberleft;
1430    double threshold = n*(n-1.0)*0.5*goDense_;
1431    if ( left >= threshold) 
1432      break;
1433    numberleft--;
1434  }
1435  //iRow=numberRows_;
1436  int nDense = numberRows_-iRow;
1437#define DENSE_THRESHOLD 8
1438  // don't do if dense columns
1439  if (nDense>=DENSE_THRESHOLD&&!dense_) {
1440    printf("Going dense for last %d rows\n",nDense);
1441    // make sure we don't disturb any indices
1442    CoinBigIndex k=0;
1443    for (int jRow=0;jRow<iRow;jRow++) {
1444      int nz=choleskyStart_[jRow+1]-choleskyStart_[jRow];
1445      k=CoinMax(k,indexStart_[jRow]+nz);
1446    }
1447    indexStart_[iRow]=k;
1448    int j;
1449    for (j=iRow+1;j<numberRows_;j++) {
1450      choleskyRow_[k++]=j;
1451      indexStart_[j]=k;
1452    }
1453    sizeIndex_=k;
1454    assert (k<=sizeFactor_); // can't happen with any reasonable defaults
1455    k=choleskyStart_[iRow];
1456    for (j=iRow+1;j<=numberRows_;j++) {
1457      k += numberRows_-j;
1458      choleskyStart_[j]=k;
1459    }
1460    // allow for blocked dense
1461    ClpCholeskyDense dense;
1462    sizeFactor_=choleskyStart_[iRow]+dense.space(nDense);
1463    firstDense_=iRow;
1464    if (doKKT_) {
1465      // redo permute so negative ones first
1466      int putN=firstDense_;
1467      int putP=0;
1468      int numberRowsModel = model_->numberRows();
1469      int numberColumns = model_->numberColumns();
1470      int numberTotal = numberColumns + numberRowsModel;
1471      for (iRow=firstDense_;iRow<numberRows_;iRow++) {
1472        int originalRow=permute_[iRow];
1473        if (originalRow<numberTotal)
1474          permute_[putN++]=originalRow;
1475        else
1476          permuteInverse_[putP++]=originalRow;
1477      }
1478      for (iRow=putN;iRow<numberRows_;iRow++) {
1479        permute_[iRow]=permuteInverse_[iRow-putN];
1480      }
1481      for (iRow=0;iRow<numberRows_;iRow++) {
1482        permuteInverse_[permute_[iRow]]=iRow;
1483      }
1484    }
1485  }
1486  // Clean up clique info
1487  for (iRow=0;iRow<numberRows_;iRow++)
1488    clique_[iRow]=0;
1489  int lastClique=-1;
1490  bool inClique=false;
1491  for (iRow=1;iRow<firstDense_;iRow++) {
1492    int sizeLast = choleskyStart_[iRow]-choleskyStart_[iRow-1];
1493    int sizeThis = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1494    if (indexStart_[iRow]==indexStart_[iRow-1]+1&&
1495        sizeThis==sizeLast-1&&
1496        sizeThis) {
1497      // in clique
1498      if (!inClique) {
1499        inClique=true;
1500        lastClique=iRow-1;
1501      }
1502    } else if (inClique) {
1503      int sizeClique=iRow-lastClique;
1504      for (int i=lastClique;i<iRow;i++) {
1505        clique_[i]=sizeClique;
1506        sizeClique--;
1507      }
1508      inClique=false;
1509    }
1510  }
1511  if (inClique) {
1512    int sizeClique=iRow-lastClique;
1513    for (int i=lastClique;i<iRow;i++) {
1514      clique_[i]=sizeClique;
1515      sizeClique--;
1516    }
1517  }
1518  //for (iRow=0;iRow<numberRows_;iRow++)
1519  //clique_[iRow]=0;
1520}
1521/* Factorize - filling in rowsDropped and returning number dropped */
1522int 
1523ClpCholeskyBase::factorize(const CoinWorkDouble * diagonal, int * rowsDropped) 
1524{
1525  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
1526  const int * columnLength = model_->clpMatrix()->getVectorLengths();
1527  const int * row = model_->clpMatrix()->getIndices();
1528  const double * element = model_->clpMatrix()->getElements();
1529  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
1530  const int * rowLength = rowCopy_->getVectorLengths();
1531  const int * column = rowCopy_->getIndices();
1532  const double * elementByRow = rowCopy_->getElements();
1533  int numberColumns=model_->clpMatrix()->getNumCols();
1534  //perturbation
1535  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
1536  //perturbation=perturbation*perturbation*100000000.0;
1537  if (perturbation>1.0) {
1538#ifdef COIN_DEVELOP
1539    //if (model_->model()->logLevel()&4)
1540      std::cout <<"large perturbation "<<perturbation<<std::endl;
1541#endif
1542    perturbation=CoinSqrt(perturbation);
1543    perturbation=1.0;
1544  }
1545  int iRow;
1546  int iColumn;
1547  longDouble * work = workDouble_;
1548  CoinZeroN(work,numberRows_);
1549  int newDropped=0;
1550  CoinWorkDouble largest=1.0;
1551  CoinWorkDouble smallest=COIN_DBL_MAX;
1552  int numberDense=0;
1553  if (!doKKT_) {
1554    const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
1555    if (dense_)
1556      numberDense=dense_->numberRows();
1557    if (whichDense_) {
1558      longDouble * denseDiagonal = dense_->diagonal();
1559      longDouble * dense = denseColumn_;
1560      int iDense=0;
1561      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1562        if (whichDense_[iColumn]) {
1563          CoinZeroN(dense,numberRows_);
1564          CoinBigIndex start=columnStart[iColumn];
1565          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1566          if (diagonal[iColumn]) {
1567            denseDiagonal[iDense++]=1.0/diagonal[iColumn];
1568            for (CoinBigIndex j=start;j<end;j++) {
1569              int jRow=row[j];
1570              dense[jRow] = element[j];
1571            }
1572          } else {
1573            denseDiagonal[iDense++]=1.0;
1574          }
1575          dense += numberRows_;
1576        }
1577      }
1578    }
1579    CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
1580    delta2 *= delta2;
1581    // largest in initial matrix
1582    CoinWorkDouble largest2=1.0e-20;
1583    for (iRow=0;iRow<numberRows_;iRow++) {
1584      longDouble * put = sparseFactor_+choleskyStart_[iRow];
1585      int * which = choleskyRow_+indexStart_[iRow];
1586      int iOriginalRow = permute_[iRow];
1587      int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1588      if (!rowLength[iOriginalRow])
1589        rowsDropped_[iOriginalRow]=1;
1590      if (!rowsDropped_[iOriginalRow]) {
1591        CoinBigIndex startRow=rowStart[iOriginalRow];
1592        CoinBigIndex endRow=rowStart[iOriginalRow]+rowLength[iOriginalRow];
1593        work[iRow] = diagonalSlack[iOriginalRow]+delta2;
1594        for (CoinBigIndex k=startRow;k<endRow;k++) {
1595          int iColumn=column[k];
1596          if (!whichDense_||!whichDense_[iColumn]) {
1597            CoinBigIndex start=columnStart[iColumn];
1598            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1599            CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
1600            for (CoinBigIndex j=start;j<end;j++) {
1601              int jRow=row[j];
1602              int jNewRow = permuteInverse_[jRow];
1603              if (jNewRow>=iRow&&!rowsDropped_[jRow]) {
1604                CoinWorkDouble value=element[j]*multiplier;
1605                work[jNewRow] += value;
1606              }
1607            }
1608          }
1609        }
1610        diagonal_[iRow]=work[iRow];
1611        largest2 = CoinMax(largest2,CoinAbs(work[iRow]));
1612        work[iRow]=0.0;
1613        int j;
1614        for (j=0;j<number;j++) {
1615          int jRow = which[j];
1616          put[j]=work[jRow];
1617          largest2 = CoinMax(largest2,CoinAbs(work[jRow]));
1618          work[jRow]=0.0;
1619        }
1620      } else {
1621        // dropped
1622        diagonal_[iRow]=1.0;
1623        int j;
1624        for (j=1;j<number;j++) {
1625          put[j]=0.0;
1626        }
1627      }
1628    }
1629    //check sizes
1630    largest2*=1.0e-20;
1631    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
1632    int numberDroppedBefore=0;
1633    for (iRow=0;iRow<numberRows_;iRow++) {
1634      int dropped=rowsDropped_[iRow];
1635      // Move to int array
1636      rowsDropped[iRow]=dropped;
1637      if (!dropped) {
1638        CoinWorkDouble diagonal = diagonal_[iRow];
1639        if (diagonal>largest2) {
1640          diagonal_[iRow]=diagonal+perturbation;
1641        } else {
1642          diagonal_[iRow]=diagonal+perturbation;
1643          rowsDropped[iRow]=2;
1644          numberDroppedBefore++;
1645          //printf("dropped - small diagonal %g\n",diagonal);
1646        } 
1647      } 
1648    }
1649    doubleParameters_[10]=CoinMax(1.0e-20,largest);
1650    integerParameters_[20]=0;
1651    doubleParameters_[3]=0.0;
1652    doubleParameters_[4]=COIN_DBL_MAX;
1653    integerParameters_[34]=0; // say all must be positive
1654    factorizePart2(rowsDropped);
1655    newDropped=integerParameters_[20]+numberDroppedBefore;
1656    largest=doubleParameters_[3];
1657    smallest=doubleParameters_[4];
1658    if (model_->messageHandler()->logLevel()>1) 
1659      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
1660    choleskyCondition_=largest/smallest;
1661    if (whichDense_) {
1662      int i;
1663      for ( i=0;i<numberRows_;i++) {
1664        assert (diagonal_[i]>=0.0);
1665        diagonal_[i]=CoinSqrt(diagonal_[i]);
1666      }
1667      // Update dense columns (just L)
1668      // Zero out dropped rows
1669      for (i=0;i<numberDense;i++) {
1670        longDouble * a = denseColumn_+i*numberRows_;
1671        for (int j=0;j<numberRows_;j++) {
1672          if (rowsDropped[j])
1673            a[j]=0.0;
1674        }
1675        for (i=0;i<numberRows_;i++) {
1676          int iRow = permute_[i];
1677          workDouble_[i] = a[iRow];
1678        }
1679        for (i=0;i<numberRows_;i++) {
1680          CoinWorkDouble value=workDouble_[i];
1681          CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
1682          CoinBigIndex j;
1683          for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
1684            int iRow = choleskyRow_[j+offset];
1685            workDouble_[iRow] -= sparseFactor_[j]*value;
1686          }
1687        }
1688        for (i=0;i<numberRows_;i++) {
1689          int iRow = permute_[i];
1690          a[iRow]=workDouble_[i]*diagonal_[i];
1691        }
1692      }
1693      dense_->resetRowsDropped();
1694      longDouble * denseBlob = dense_->aMatrix();
1695      longDouble * denseDiagonal = dense_->diagonal();
1696      // Update dense matrix
1697      for (i=0;i<numberDense;i++) {
1698        const longDouble * a = denseColumn_+i*numberRows_;
1699        // do diagonal
1700        CoinWorkDouble value = denseDiagonal[i];
1701        const longDouble * b = denseColumn_+i*numberRows_;
1702        for (int k=0;k<numberRows_;k++) 
1703          value += a[k]*b[k];
1704        denseDiagonal[i]=value;
1705        for (int j=i+1;j<numberDense;j++) {
1706          CoinWorkDouble value = 0.0;
1707          const longDouble * b = denseColumn_+j*numberRows_;
1708          for (int k=0;k<numberRows_;k++) 
1709            value += a[k]*b[k];
1710          *denseBlob=value;
1711          denseBlob++;
1712        }
1713      }
1714      // dense cholesky (? long double)
1715      int * dropped = new int [numberDense];
1716      dense_->factorizePart2(dropped);
1717      delete [] dropped;
1718    }
1719    // try allowing all every time
1720    //printf("trying ?\n");
1721    //for (iRow=0;iRow<numberRows_;iRow++) {
1722    //rowsDropped[iRow]=0;
1723    //rowsDropped_[iRow]=0;
1724    //}
1725    bool cleanCholesky;
1726    //if (model_->numberIterations()<20||(model_->numberIterations()&1)==0)
1727    if (model_->numberIterations()<2000) 
1728      cleanCholesky=true;
1729    else 
1730      cleanCholesky=false;
1731    if (cleanCholesky) {
1732      //drop fresh makes some formADAT easier
1733      if (newDropped||numberRowsDropped_) {
1734        newDropped=0;
1735        for (int i=0;i<numberRows_;i++) {
1736          char dropped = static_cast<char>(rowsDropped[i]);
1737          rowsDropped_[i]=dropped;
1738          rowsDropped_[i]=0;
1739          if (dropped==2) {
1740            //dropped this time
1741            rowsDropped[newDropped++]=i;
1742            rowsDropped_[i]=0;
1743          } 
1744        } 
1745        numberRowsDropped_=newDropped;
1746        newDropped=-(2+newDropped);
1747      } 
1748    } else {
1749      if (newDropped) {
1750        newDropped=0;
1751        for (int i=0;i<numberRows_;i++) {
1752          char dropped = static_cast<char>(rowsDropped[i]);
1753          rowsDropped_[i]=dropped;
1754          if (dropped==2) {
1755            //dropped this time
1756            rowsDropped[newDropped++]=i;
1757            rowsDropped_[i]=1;
1758          } 
1759        } 
1760      } 
1761      numberRowsDropped_+=newDropped;
1762      if (numberRowsDropped_&&0) {
1763        std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
1764          numberRowsDropped_<<" dropped)";
1765        if (newDropped) {
1766          std::cout<<" ( "<<newDropped<<" dropped this time)";
1767        } 
1768        std::cout<<std::endl;
1769      } 
1770    }
1771  } else {
1772    //KKT
1773    CoinPackedMatrix * quadratic = NULL;
1774    ClpQuadraticObjective * quadraticObj = 
1775      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
1776    if (quadraticObj) 
1777      quadratic = quadraticObj->quadraticObjective();
1778    // matrix
1779    int numberRowsModel = model_->numberRows();
1780    int numberColumns = model_->numberColumns();
1781    int numberTotal = numberColumns + numberRowsModel;
1782    // temp
1783    bool permuted=false;
1784    for (iRow=0;iRow<numberRows_;iRow++) {
1785      if (permute_[iRow]!=iRow) {
1786        permuted=true;
1787        break;
1788      }
1789    }
1790    // but fake it
1791    for (iRow=0;iRow<numberRows_;iRow++) {
1792      //permute_[iRow]=iRow; // force no permute
1793      //permuteInverse_[permute_[iRow]]=iRow;
1794    }
1795    if (permuted) {
1796      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
1797      delta2 *= delta2;
1798      // Need to permute - ugly
1799      if (!quadratic) {
1800        for (iRow=0;iRow<numberRows_;iRow++) {
1801          longDouble * put = sparseFactor_+choleskyStart_[iRow];
1802          int * which = choleskyRow_+indexStart_[iRow];
1803          int iOriginalRow = permute_[iRow];
1804          if (iOriginalRow<numberColumns) {
1805            iColumn=iOriginalRow;
1806            CoinWorkDouble value = diagonal[iColumn];
1807            if (CoinAbs(value)>1.0e-100) {
1808              value = 1.0/value;
1809              largest = CoinMax(largest,CoinAbs(value));
1810              diagonal_[iRow] = -value;
1811              CoinBigIndex start=columnStart[iColumn];
1812              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1813              for (CoinBigIndex j=start;j<end;j++) {
1814                int kRow = row[j]+numberTotal;
1815                kRow=permuteInverse_[kRow];
1816                if (kRow>iRow) {
1817                  work[kRow]=element[j];
1818                  largest = CoinMax(largest,CoinAbs(element[j]));
1819                }
1820              }
1821            } else {
1822              diagonal_[iRow] = -value;
1823            }
1824          } else if (iOriginalRow<numberTotal) {
1825            CoinWorkDouble value = diagonal[iOriginalRow];
1826            if (CoinAbs(value)>1.0e-100) {
1827              value = 1.0/value;
1828              largest = CoinMax(largest,CoinAbs(value));
1829            } else {
1830              value = 1.0e100;
1831            }
1832            diagonal_[iRow] = -value;
1833            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1834            if (kRow>iRow) 
1835              work[kRow]=-1.0;
1836          } else {
1837            diagonal_[iRow]=delta2;
1838            int kRow = iOriginalRow-numberTotal;
1839            CoinBigIndex start=rowStart[kRow];
1840            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1841            for (CoinBigIndex j=start;j<end;j++) {
1842              int jRow = column[j];
1843              int jNewRow = permuteInverse_[jRow];
1844              if (jNewRow>iRow) {
1845                work[jNewRow]=elementByRow[j];
1846                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
1847              }
1848            }
1849            // slack - should it be permute
1850            kRow = permuteInverse_[kRow+numberColumns];
1851            if (kRow>iRow)
1852              work[kRow]=-1.0;
1853          }
1854          CoinBigIndex j;
1855          int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1856          for (j=0;j<number;j++) {
1857            int jRow = which[j];
1858            put[j]=work[jRow];
1859            work[jRow]=0.0;
1860          }
1861        }
1862      } else {
1863        // quadratic
1864        const int * columnQuadratic = quadratic->getIndices();
1865        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1866        const int * columnQuadraticLength = quadratic->getVectorLengths();
1867        const double * quadraticElement = quadratic->getElements();
1868        for (iRow=0;iRow<numberRows_;iRow++) {
1869          longDouble * put = sparseFactor_+choleskyStart_[iRow];
1870          int * which = choleskyRow_+indexStart_[iRow];
1871          int iOriginalRow = permute_[iRow];
1872          if (iOriginalRow<numberColumns) {
1873            CoinBigIndex j;
1874            iColumn=iOriginalRow;
1875            CoinWorkDouble value = diagonal[iColumn];
1876            if (CoinAbs(value)>1.0e-100) {
1877              value = 1.0/value;
1878              for (j=columnQuadraticStart[iColumn];
1879                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1880                int jColumn = columnQuadratic[j];
1881                int jNewColumn = permuteInverse_[jColumn];
1882                if (jNewColumn>iRow) {
1883                  work[jNewColumn]=-quadraticElement[j];
1884                } else if (iColumn==jColumn) {
1885                  value += quadraticElement[j];
1886                }
1887              }
1888              largest = CoinMax(largest,CoinAbs(value));
1889              diagonal_[iRow] = -value;
1890              CoinBigIndex start=columnStart[iColumn];
1891              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1892              for (j=start;j<end;j++) {
1893                int kRow = row[j]+numberTotal;
1894                kRow=permuteInverse_[kRow];
1895                if (kRow>iRow) {
1896                  work[kRow]=element[j];
1897                  largest = CoinMax(largest,CoinAbs(element[j]));
1898                }
1899              }
1900            } else {
1901              diagonal_[iRow] = -value;
1902            }
1903          } else if (iOriginalRow<numberTotal) {
1904            CoinWorkDouble value = diagonal[iOriginalRow];
1905            if (CoinAbs(value)>1.0e-100) {
1906              value = 1.0/value;
1907              largest = CoinMax(largest,CoinAbs(value));
1908            } else {
1909              value = 1.0e100;
1910            }
1911            diagonal_[iRow] = -value;
1912            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1913            if (kRow>iRow) 
1914              work[kRow]=-1.0;
1915          } else {
1916            diagonal_[iRow]=delta2;
1917            int kRow = iOriginalRow-numberTotal;
1918            CoinBigIndex start=rowStart[kRow];
1919            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1920            for (CoinBigIndex j=start;j<end;j++) {
1921              int jRow = column[j];
1922              int jNewRow = permuteInverse_[jRow];
1923              if (jNewRow>iRow) {
1924                work[jNewRow]=elementByRow[j];
1925                largest = CoinMax(largest,CoinAbs(elementByRow[j]));
1926              }
1927            }
1928            // slack - should it be permute
1929            kRow = permuteInverse_[kRow+numberColumns];
1930            if (kRow>iRow)
1931              work[kRow]=-1.0;
1932          }
1933          CoinBigIndex j;
1934          int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1935          for (j=0;j<number;j++) {
1936            int jRow = which[j];
1937            put[j]=work[jRow];
1938            work[jRow]=0.0;
1939          }
1940          for (j=0;j<numberRows_;j++)
1941            assert (!work[j]);
1942        }
1943      }
1944    } else {
1945      if (!quadratic) {
1946        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1947          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
1948          int * which = choleskyRow_+indexStart_[iColumn];
1949          CoinWorkDouble value = diagonal[iColumn];
1950          if (CoinAbs(value)>1.0e-100) {
1951            value = 1.0/value;
1952            largest = CoinMax(largest,CoinAbs(value));
1953            diagonal_[iColumn] = -value;
1954            CoinBigIndex start=columnStart[iColumn];
1955            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1956            for (CoinBigIndex j=start;j<end;j++) {
1957              //choleskyRow_[numberElements]=row[j]+numberTotal;
1958              //sparseFactor_[numberElements++]=element[j];
1959              work[row[j]+numberTotal]=element[j];
1960              largest = CoinMax(largest,CoinAbs(element[j]));
1961            }
1962          } else {
1963            diagonal_[iColumn] = -value;
1964          }
1965          CoinBigIndex j;
1966          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
1967          for (j=0;j<number;j++) {
1968            int jRow = which[j];
1969            put[j]=work[jRow];
1970            work[jRow]=0.0;
1971          }
1972        }
1973      } else {
1974        // Quadratic
1975        const int * columnQuadratic = quadratic->getIndices();
1976        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1977        const int * columnQuadraticLength = quadratic->getVectorLengths();
1978        const double * quadraticElement = quadratic->getElements();
1979        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1980          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
1981          int * which = choleskyRow_+indexStart_[iColumn];
1982          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
1983          CoinWorkDouble value = diagonal[iColumn];
1984          CoinBigIndex j;
1985          if (CoinAbs(value)>1.0e-100) {
1986            value = 1.0/value;
1987            for (j=columnQuadraticStart[iColumn];
1988                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1989              int jColumn = columnQuadratic[j];
1990              if (jColumn>iColumn) {
1991                work[jColumn]=-quadraticElement[j];
1992              } else if (iColumn==jColumn) {
1993                value += quadraticElement[j];
1994              }
1995            }
1996            largest = CoinMax(largest,CoinAbs(value));
1997            diagonal_[iColumn] = -value;
1998            CoinBigIndex start=columnStart[iColumn];
1999            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
2000            for (j=start;j<end;j++) {
2001              work[row[j]+numberTotal]=element[j];
2002              largest = CoinMax(largest,CoinAbs(element[j]));
2003            }
2004            for (j=0;j<number;j++) {
2005              int jRow = which[j];
2006              put[j]=work[jRow];
2007              work[jRow]=0.0;
2008            }
2009          } else {
2010            value = 1.0e100;
2011            diagonal_[iColumn] = -value;
2012            for (j=0;j<number;j++) {
2013              int jRow = which[j];
2014              put[j]=work[jRow];
2015            }
2016          }
2017        }
2018      }
2019      // slacks
2020      for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
2021        longDouble * put = sparseFactor_+choleskyStart_[iColumn];
2022        int * which = choleskyRow_+indexStart_[iColumn];
2023        CoinWorkDouble value = diagonal[iColumn];
2024        if (CoinAbs(value)>1.0e-100) {
2025          value = 1.0/value;
2026          largest = CoinMax(largest,CoinAbs(value));
2027        } else {
2028          value = 1.0e100;
2029        }
2030        diagonal_[iColumn] = -value;
2031        work[iColumn-numberColumns+numberTotal]=-1.0;
2032        CoinBigIndex j;
2033        int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
2034        for (j=0;j<number;j++) {
2035          int jRow = which[j];
2036          put[j]=work[jRow];
2037          work[jRow]=0.0;
2038        }
2039      }
2040      // Finish diagonal
2041      CoinWorkDouble delta2 = model_->delta(); // add delta*delta to bottom
2042      delta2 *= delta2;
2043      for (iRow=0;iRow<numberRowsModel;iRow++) {
2044        longDouble * put = sparseFactor_+choleskyStart_[iRow+numberTotal];
2045        diagonal_[iRow+numberTotal]=delta2;
2046        CoinBigIndex j;
2047        int number = choleskyStart_[iRow+numberTotal+1]-choleskyStart_[iRow+numberTotal];
2048        for (j=0;j<number;j++) {
2049          put[j]=0.0;
2050        }
2051      }
2052    }
2053    //check sizes
2054    largest*=1.0e-20;
2055    largest = CoinMin(largest,CHOL_SMALL_VALUE);
2056    doubleParameters_[10]=CoinMax(1.0e-20,largest);
2057    integerParameters_[20]=0;
2058    doubleParameters_[3]=0.0;
2059    doubleParameters_[4]=COIN_DBL_MAX;
2060    // Set up LDL cutoff
2061    integerParameters_[34]=numberTotal;
2062    // KKT
2063    int * rowsDropped2 = new int[numberRows_];
2064    CoinZeroN(rowsDropped2,numberRows_);
2065    factorizePart2(rowsDropped2);
2066    newDropped=integerParameters_[20];
2067    largest=doubleParameters_[3];
2068    smallest=doubleParameters_[4];
2069    if (model_->messageHandler()->logLevel()>1) 
2070      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
2071    choleskyCondition_=largest/smallest;
2072    // Should save adjustments in ..R_
2073    int n1=0,n2=0;
2074    CoinWorkDouble * primalR = model_->primalR();
2075    CoinWorkDouble * dualR = model_->dualR();
2076    for (iRow=0;iRow<numberTotal;iRow++) { 
2077      if (rowsDropped2[iRow]) {
2078        n1++;
2079        //printf("row region1 %d dropped\n",iRow);
2080        //rowsDropped_[iRow]=1;
2081        rowsDropped_[iRow]=0;
2082        primalR[iRow]=doubleParameters_[20];
2083      } else {
2084        rowsDropped_[iRow]=0;
2085        primalR[iRow]=0.0;
2086      }
2087    }
2088    for (;iRow<numberRows_;iRow++) {
2089      if (rowsDropped2[iRow]) {
2090        n2++;
2091        //printf("row region2 %d dropped\n",iRow);
2092        //rowsDropped_[iRow]=1;
2093        rowsDropped_[iRow]=0;
2094        dualR[iRow-numberTotal]=doubleParameters_[34];
2095      } else {
2096        rowsDropped_[iRow]=0;
2097        dualR[iRow-numberTotal]=0.0;
2098      }
2099    }
2100    delete [] rowsDropped2;
2101  }
2102  status_=0;
2103  return newDropped;
2104}
2105/* Factorize - filling in rowsDropped and returning number dropped
2106   in integerParam.
2107*/
2108void 
2109ClpCholeskyBase::factorizePart2(int * rowsDropped) 
2110{
2111  CoinWorkDouble largest=doubleParameters_[3];
2112  CoinWorkDouble smallest=doubleParameters_[4];
2113  int numberDropped=integerParameters_[20];
2114  // probably done before
2115  largest=0.0;
2116  smallest=COIN_DBL_MAX;
2117  numberDropped=0;
2118  double dropValue = doubleParameters_[10];
2119  int firstPositive=integerParameters_[34];
2120  longDouble * d = ClpCopyOfArray(diagonal_,numberRows_);
2121  int iRow;
2122  // minimum size before clique done
2123  //#define MINCLIQUE INT_MAX
2124#define MINCLIQUE 3
2125  longDouble * work = workDouble_;
2126  CoinBigIndex * first = workInteger_;
2127 
2128  for (iRow=0;iRow<numberRows_;iRow++) {
2129    link_[iRow]=-1;
2130    work[iRow]=0.0;
2131    first[iRow]=choleskyStart_[iRow];
2132  }
2133
2134  int lastClique=-1;
2135  bool inClique=false;
2136  bool newClique=false;
2137  bool endClique=false;
2138  int lastRow=0;
2139  int cliqueSize=0;
2140  CoinBigIndex cliquePointer=0;
2141  int nextRow2=-1;
2142 
2143  for (iRow=0;iRow<firstDense_+1;iRow++) {
2144    if (iRow<firstDense_) {
2145      endClique=false;
2146      if (clique_[iRow]>0) {
2147        // this is in a clique
2148        inClique=true;
2149        if (clique_[iRow]>lastClique) {
2150          // new Clique
2151          newClique=true;
2152          // If we have clique going then signal to do old one
2153          endClique=(lastClique>0);
2154        } else {
2155          // Still in clique
2156          newClique=false;
2157        }
2158      } else {
2159        // not in clique
2160        inClique=false;
2161        newClique=false;
2162        // If we have clique going then signal to do old one
2163        endClique=(lastClique>0);
2164      }
2165      lastClique=clique_[iRow];
2166    } else if (inClique) {
2167      // Finish off
2168      endClique=true;
2169    } else {
2170      break;
2171    }
2172    if (endClique) {
2173      // We have just finished updating a clique - do block pivot and clean up
2174      int jRow;
2175      for ( jRow=lastRow;jRow<iRow;jRow++) {
2176        int jCount = jRow-lastRow;
2177        CoinWorkDouble diagonalValue = diagonal_[jRow];
2178        CoinBigIndex start=choleskyStart_[jRow];
2179        CoinBigIndex end=choleskyStart_[jRow+1];
2180        for (int kRow=lastRow;kRow<jRow;kRow++) {
2181          jCount--;
2182          CoinBigIndex get = choleskyStart_[kRow]+jCount;
2183          CoinWorkDouble a_jk = sparseFactor_[get];
2184          CoinWorkDouble value1 = d[kRow]*a_jk;
2185          diagonalValue -= a_jk*value1;
2186          for (CoinBigIndex j=start;j<end;j++)
2187            sparseFactor_[j] -= value1*sparseFactor_[++get];
2188        }
2189        // check
2190        int originalRow = permute_[jRow];
2191        if (originalRow<firstPositive) {
2192          // must be negative
2193          if (diagonalValue<=-dropValue) {
2194            smallest = CoinMin(smallest,-diagonalValue);
2195            largest = CoinMax(largest,-diagonalValue);
2196            d[jRow]=diagonalValue;
2197            diagonalValue = 1.0/diagonalValue;
2198          } else {
2199            rowsDropped[originalRow]=2;
2200            d[jRow]=-1.0e100;
2201            diagonalValue=0.0;
2202            integerParameters_[20]++;
2203          }
2204        } else {
2205          // must be positive
2206          if (diagonalValue>=dropValue) {
2207            smallest = CoinMin(smallest,diagonalValue);
2208            largest = CoinMax(largest,diagonalValue);
2209            d[jRow]=diagonalValue;
2210            diagonalValue = 1.0/diagonalValue;
2211          } else {
2212            rowsDropped[originalRow]=2;
2213            d[jRow]=1.0e100;
2214            diagonalValue=0.0;
2215            integerParameters_[20]++;
2216          }
2217        }
2218        diagonal_[jRow]=diagonalValue;
2219        for (CoinBigIndex j=start;j<end;j++) {
2220          sparseFactor_[j] *= diagonalValue;
2221        }
2222      }
2223      if (nextRow2>=0) {
2224        for ( jRow=lastRow;jRow<iRow-1;jRow++) {
2225          link_[jRow]=jRow+1;
2226        }
2227        link_[iRow-1]=link_[nextRow2];
2228        link_[nextRow2]=lastRow;
2229      }
2230    }
2231    if (iRow==firstDense_)
2232      break; // we were just cleaning up
2233    if (newClique) {
2234      // initialize new clique
2235      lastRow=iRow;
2236      cliquePointer=choleskyStart_[iRow];
2237      cliqueSize=choleskyStart_[iRow+1]-cliquePointer+1;
2238    }
2239    // for each column L[*,kRow] that affects L[*,iRow]
2240    CoinWorkDouble diagonalValue=diagonal_[iRow];
2241    int nextRow = link_[iRow];
2242    int kRow=0;
2243    while (1) {
2244      kRow=nextRow;
2245      if (kRow<0)
2246        break; // out of loop
2247      nextRow=link_[kRow];
2248      // Modify by outer product of L[*,irow] by L[*,krow] from first
2249      CoinBigIndex k=first[kRow];
2250      CoinBigIndex end=choleskyStart_[kRow+1];
2251      assert(k<end);
2252      CoinWorkDouble a_ik=sparseFactor_[k++];
2253      CoinWorkDouble value1=d[kRow]*a_ik;
2254      // update first
2255      first[kRow]=k;
2256      diagonalValue -= value1*a_ik;
2257      CoinBigIndex offset = indexStart_[kRow]-choleskyStart_[kRow];
2258      if (k<end) {
2259        int jRow = choleskyRow_[k+offset];
2260        if (clique_[kRow]<MINCLIQUE) {
2261          link_[kRow]=link_[jRow];
2262          link_[jRow]=kRow;
2263          for (;k<end;k++) {
2264            int jRow = choleskyRow_[k+offset];
2265            work[jRow] += sparseFactor_[k]*value1;
2266          }
2267        } else {
2268          // Clique
2269          CoinBigIndex currentIndex = k+offset;
2270          int linkSave=link_[jRow];
2271          link_[jRow]=kRow;
2272          work[kRow]=value1; // ? or a_jk
2273          int last = kRow+clique_[kRow];
2274          for (int kkRow=kRow+1;kkRow<last;kkRow++) {
2275            CoinBigIndex j=first[kkRow];
2276            //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
2277            CoinWorkDouble a = sparseFactor_[j];
2278            CoinWorkDouble dValue = d[kkRow]*a;
2279            diagonalValue -= a*dValue;
2280            work[kkRow]=dValue;
2281            first[kkRow]++;
2282            link_[kkRow-1]=kkRow;
2283          }
2284          nextRow = link_[last-1];
2285          link_[last-1]=linkSave;
2286          int length = end-k;
2287          for (int i=0;i<length;i++) {
2288            int lRow = choleskyRow_[currentIndex++];
2289            CoinWorkDouble t0 = work[lRow];
2290            for (int kkRow=kRow;kkRow<last;kkRow++) {
2291              CoinBigIndex j = first[kkRow]+i;
2292              t0 += work[kkRow]*sparseFactor_[j];
2293            }
2294            work[lRow]=t0;
2295          }
2296        }
2297      }
2298    }
2299    // Now apply
2300    if (inClique) {
2301      // in clique
2302      diagonal_[iRow]=diagonalValue;
2303      CoinBigIndex start=choleskyStart_[iRow];
2304      CoinBigIndex end=choleskyStart_[iRow+1];
2305      CoinBigIndex currentIndex = indexStart_[iRow];
2306      nextRow2=-1;
2307      CoinBigIndex get=start+clique_[iRow]-1;
2308      if (get<end) {
2309        nextRow2 = choleskyRow_[currentIndex+get-start];
2310        first[iRow]=get;
2311      }
2312      for (CoinBigIndex j=start;j<end;j++) {
2313        int kRow = choleskyRow_[currentIndex++];
2314        sparseFactor_[j] -= work[kRow]; // times?
2315        work[kRow]=0.0;
2316      }
2317    } else {
2318      // not in clique
2319      int originalRow = permute_[iRow];
2320      if (originalRow<firstPositive) {
2321        // must be negative
2322        if (diagonalValue<=-dropValue) {
2323          smallest = CoinMin(smallest,-diagonalValue);
2324          largest = CoinMax(largest,-diagonalValue);
2325          d[iRow]=diagonalValue;
2326          diagonalValue = 1.0/diagonalValue;
2327        } else {
2328          rowsDropped[originalRow]=2;
2329          d[iRow]=-1.0e100;
2330          diagonalValue=0.0;
2331          integerParameters_[20]++;
2332        }
2333      } else {
2334        // must be positive
2335        if (diagonalValue>=dropValue) {
2336          smallest = CoinMin(smallest,diagonalValue);
2337          largest = CoinMax(largest,diagonalValue);
2338          d[iRow]=diagonalValue;
2339          diagonalValue = 1.0/diagonalValue;
2340        } else {
2341          rowsDropped[originalRow]=2;
2342          d[iRow]=1.0e100;
2343          diagonalValue=0.0;
2344          integerParameters_[20]++;
2345        }
2346      }
2347      diagonal_[iRow]=diagonalValue;
2348      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
2349      CoinBigIndex start = choleskyStart_[iRow];
2350      CoinBigIndex end = choleskyStart_[iRow+1];
2351      assert (first[iRow]==start);
2352      if (start<end) {
2353        int nextRow = choleskyRow_[start+offset];
2354        link_[iRow]=link_[nextRow];
2355        link_[nextRow]=iRow;
2356        for (CoinBigIndex j=start;j<end;j++) {
2357          int jRow = choleskyRow_[j+offset];
2358          CoinWorkDouble value = sparseFactor_[j]-work[jRow];
2359          work[jRow]=0.0;
2360          sparseFactor_[j]= diagonalValue*value;
2361        }
2362      }
2363    }
2364  }
2365  if (firstDense_<numberRows_) {
2366    // do dense
2367    // update dense part
2368    updateDense(d,work,first);
2369    ClpCholeskyDense dense;
2370    // just borrow space
2371    int nDense = numberRows_-firstDense_;
2372    if (doKKT_) {
2373      for (iRow=firstDense_;iRow<numberRows_;iRow++) {
2374        int originalRow=permute_[iRow];
2375        if (originalRow>=firstPositive) {
2376          firstPositive=iRow-firstDense_;
2377          break;
2378        }
2379      }
2380    }
2381    dense.reserveSpace(this,nDense);
2382    int * dropped = new int[nDense];
2383    memset(dropped,0,nDense*sizeof(int));
2384    dense.setDoubleParameter(3,largest);
2385    dense.setDoubleParameter(4,smallest);
2386    dense.setDoubleParameter(10,dropValue);
2387    dense.setIntegerParameter(20,0);
2388    dense.setIntegerParameter(34,firstPositive);
2389    dense.factorizePart2(dropped);
2390    largest=dense.getDoubleParameter(3);
2391    smallest=dense.getDoubleParameter(4);
2392    integerParameters_[20]+=dense.getIntegerParameter(20);
2393    for (iRow=firstDense_;iRow<numberRows_;iRow++) {
2394      int originalRow=permute_[iRow];
2395      rowsDropped[originalRow]=dropped[iRow-firstDense_];
2396    }
2397    delete [] dropped;
2398  }
2399  delete [] d;
2400  doubleParameters_[3]=largest;
2401  doubleParameters_[4]=smallest;
2402  return;
2403}
2404// Updates dense part (broken out for profiling)
2405void ClpCholeskyBase::updateDense(longDouble * d, longDouble * work, int * first)
2406{
2407  for (int iRow=0;iRow<firstDense_;iRow++) {
2408    CoinBigIndex start=first[iRow];
2409    CoinBigIndex end=choleskyStart_[iRow+1];
2410    if (start<end) {
2411      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
2412      if (clique_[iRow]<2) {
2413        CoinWorkDouble dValue=d[iRow];
2414        for (CoinBigIndex k=start;k<end;k++) {
2415          int kRow = choleskyRow_[k+offset]; 
2416          assert(kRow>=firstDense_);
2417          CoinWorkDouble a_ik=sparseFactor_[k];
2418          CoinWorkDouble value1=dValue*a_ik;
2419          diagonal_[kRow] -= value1*a_ik;
2420          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2421          for (CoinBigIndex j=k+1;j<end;j++) {
2422            int jRow=choleskyRow_[j+offset];
2423            CoinWorkDouble a_jk = sparseFactor_[j];
2424            sparseFactor_[base+jRow] -= a_jk*value1;
2425          }
2426        }
2427      } else if (clique_[iRow]<3) {
2428        // do as pair
2429        CoinWorkDouble dValue0=d[iRow];
2430        CoinWorkDouble dValue1=d[iRow+1];
2431        int offset1 = first[iRow+1]-start;
2432        // skip row
2433        iRow++;
2434        for (CoinBigIndex k=start;k<end;k++) {
2435          int kRow = choleskyRow_[k+offset]; 
2436          assert(kRow>=firstDense_);
2437          CoinWorkDouble a_ik0=sparseFactor_[k];
2438          CoinWorkDouble value0=dValue0*a_ik0;
2439          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
2440          CoinWorkDouble value1=dValue1*a_ik1;
2441          diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1;
2442          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2443          for (CoinBigIndex j=k+1;j<end;j++) {
2444            int jRow=choleskyRow_[j+offset];
2445            CoinWorkDouble a_jk0 = sparseFactor_[j];
2446            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
2447            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1;
2448          }
2449        }
2450#define MANY_REGISTERS
2451#ifdef MANY_REGISTERS
2452      } else if (clique_[iRow]==3) {
2453#else
2454      } else {
2455#endif
2456        // do as clique
2457        // maybe later get fancy on big cliques and do transpose copy
2458        // seems only worth going to 3 on Intel
2459        CoinWorkDouble dValue0=d[iRow];
2460        CoinWorkDouble dValue1=d[iRow+1];
2461        CoinWorkDouble dValue2=d[iRow+2];
2462        // get offsets and skip rows
2463        int offset1 = first[++iRow]-start;
2464        int offset2 = first[++iRow]-start;
2465        for (CoinBigIndex k=start;k<end;k++) {
2466          int kRow = choleskyRow_[k+offset]; 
2467          assert(kRow>=firstDense_);
2468          CoinWorkDouble diagonalValue=diagonal_[kRow];
2469          CoinWorkDouble a_ik0=sparseFactor_[k];
2470          CoinWorkDouble value0=dValue0*a_ik0;
2471          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
2472          CoinWorkDouble value1=dValue1*a_ik1;
2473          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
2474          CoinWorkDouble value2=dValue2*a_ik2;
2475          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2476          diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2;
2477          for (CoinBigIndex j=k+1;j<end;j++) {
2478            int jRow=choleskyRow_[j+offset];
2479            CoinWorkDouble a_jk0 = sparseFactor_[j];
2480            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
2481            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
2482            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2;
2483          }
2484        }
2485#ifdef MANY_REGISTERS
2486      } else {
2487        // do as clique
2488        // maybe later get fancy on big cliques and do transpose copy
2489        // maybe only worth going to 3 on Intel (but may have hidden registers)
2490        CoinWorkDouble dValue0=d[iRow];
2491        CoinWorkDouble dValue1=d[iRow+1];
2492        CoinWorkDouble dValue2=d[iRow+2];
2493        CoinWorkDouble dValue3=d[iRow+3];
2494        // get offsets and skip rows
2495        int offset1 = first[++iRow]-start;
2496        int offset2 = first[++iRow]-start;
2497        int offset3 = first[++iRow]-start;
2498        for (CoinBigIndex k=start;k<end;k++) {
2499          int kRow = choleskyRow_[k+offset]; 
2500          assert(kRow>=firstDense_);
2501          CoinWorkDouble diagonalValue=diagonal_[kRow];
2502          CoinWorkDouble a_ik0=sparseFactor_[k];
2503          CoinWorkDouble value0=dValue0*a_ik0;
2504          CoinWorkDouble a_ik1=sparseFactor_[k+offset1];
2505          CoinWorkDouble value1=dValue1*a_ik1;
2506          CoinWorkDouble a_ik2=sparseFactor_[k+offset2];
2507          CoinWorkDouble value2=dValue2*a_ik2;
2508          CoinWorkDouble a_ik3=sparseFactor_[k+offset3];
2509          CoinWorkDouble value3=dValue3*a_ik3;
2510          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2511          diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3);
2512          for (CoinBigIndex j=k+1;j<end;j++) {
2513            int jRow=choleskyRow_[j+offset];
2514            CoinWorkDouble a_jk0 = sparseFactor_[j];
2515            CoinWorkDouble a_jk1 = sparseFactor_[j+offset1];
2516            CoinWorkDouble a_jk2 = sparseFactor_[j+offset2];
2517            CoinWorkDouble a_jk3 = sparseFactor_[j+offset3];
2518            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3;
2519          }
2520        }
2521#endif
2522      }
2523    }
2524  }
2525}
2526/* Uses factorization to solve. */
2527void 
2528ClpCholeskyBase::solve (double * region) 
2529{
2530  if (!whichDense_) {
2531    solve(region,3);
2532  } else {
2533    // dense columns
2534    int i;
2535    solve(region,1);
2536    // do change;
2537    int numberDense = dense_->numberRows();
2538    CoinWorkDouble * change = new CoinWorkDouble[numberDense];
2539    for (i=0;i<numberDense;i++) {
2540      const longDouble * a = denseColumn_+i*numberRows_;
2541      CoinWorkDouble value =0.0;
2542      for (int iRow=0;iRow<numberRows_;iRow++) 
2543        value += a[iRow]*region[iRow];
2544      change[i]=value;
2545    }
2546    // solve
2547#if CLP_LONG_CHOLESKY>0
2548    dense_->solveLong(change);
2549#else
2550    dense_->solve(change);
2551#endif
2552    for (i=0;i<numberDense;i++) {
2553      const longDouble * a = denseColumn_+i*numberRows_;
2554      CoinWorkDouble value = change[i];
2555      for (int iRow=0;iRow<numberRows_;iRow++) 
2556        region[iRow] -= value*a[iRow];
2557    }
2558    delete [] change;
2559    // and finish off
2560    solve(region,2);
2561  }
2562}
2563/* solve - 1 just first half, 2 just second half - 3 both.
2564   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
2565*/
2566void 
2567ClpCholeskyBase::solve(double * region, int type)
2568{
2569#ifdef CLP_DEBUG
2570  double * regionX=NULL;
2571  if (sizeof(CoinWorkDouble)!=sizeof(double)&&type==3) {
2572    regionX=ClpCopyOfArray(region,numberRows_);
2573  }
2574#endif
2575  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
2576  int i;
2577  CoinBigIndex j;
2578  for (i=0;i<numberRows_;i++) {
2579    int iRow = permute_[i];
2580    work[i] = region[iRow];
2581  }
2582  switch (type) {
2583  case 1:
2584    for (i=0;i<numberRows_;i++) {
2585      CoinWorkDouble value=work[i];
2586      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2587      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2588        int iRow = choleskyRow_[j+offset];
2589        work[iRow] -= sparseFactor_[j]*value;
2590      }
2591    }
2592    for (i=0;i<numberRows_;i++) {
2593      int iRow = permute_[i];
2594      region[iRow]=work[i]*diagonal_[i];
2595    }
2596    break;
2597  case 2:
2598    for (i=numberRows_-1;i>=0;i--) {
2599      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2600      CoinWorkDouble value=work[i]*diagonal_[i];
2601      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2602        int iRow = choleskyRow_[j+offset];
2603        value -= sparseFactor_[j]*work[iRow];
2604      }
2605      work[i]=value;
2606      int iRow = permute_[i];
2607      region[iRow]=value;
2608    }
2609    break;
2610  case 3:
2611    for (i=0;i<firstDense_;i++) {
2612      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2613      CoinWorkDouble value=work[i];
2614      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2615        int iRow = choleskyRow_[j+offset];
2616        work[iRow] -= sparseFactor_[j]*value;
2617      }
2618    }
2619    if (firstDense_<numberRows_) {
2620      // do dense
2621      ClpCholeskyDense dense;
2622      // just borrow space
2623      int nDense = numberRows_-firstDense_;
2624      dense.reserveSpace(this,nDense);
2625#if CLP_LONG_CHOLESKY!=1
2626      dense.solveLong(work+firstDense_);
2627#else
2628      dense.solveLongWork(work+firstDense_);
2629#endif
2630      for (i=numberRows_-1;i>=firstDense_;i--) {
2631        CoinWorkDouble value=work[i];
2632        int iRow = permute_[i];
2633        region[iRow]=value;
2634      }
2635    }
2636    for (i=firstDense_-1;i>=0;i--) {
2637      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2638      CoinWorkDouble value=work[i]*diagonal_[i];
2639      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2640        int iRow = choleskyRow_[j+offset];
2641        value -= sparseFactor_[j]*work[iRow];
2642      }
2643      work[i]=value;
2644      int iRow = permute_[i];
2645      region[iRow]=value;
2646    }
2647    break;
2648  }
2649#ifdef CLP_DEBUG
2650  if (regionX) {
2651    longDouble * work = workDouble_;
2652    int i;
2653    CoinBigIndex j;
2654    double largestO=0.0;
2655    for (i=0;i<numberRows_;i++) {
2656      largestO = CoinMax(largestO,CoinAbs(regionX[i]));
2657    }
2658    for (i=0;i<numberRows_;i++) {
2659      int iRow = permute_[i];
2660      work[i] = regionX[iRow];
2661    }
2662    for (i=0;i<firstDense_;i++) {
2663      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2664      CoinWorkDouble value=work[i];
2665      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2666        int iRow = choleskyRow_[j+offset];
2667        work[iRow] -= sparseFactor_[j]*value;
2668      }
2669    }
2670    if (firstDense_<numberRows_) {
2671      // do dense
2672      ClpCholeskyDense dense;
2673      // just borrow space
2674      int nDense = numberRows_-firstDense_;
2675      dense.reserveSpace(this,nDense);
2676      dense.solveLong(work+firstDense_);
2677      for (i=numberRows_-1;i>=firstDense_;i--) {
2678        CoinWorkDouble value=work[i];
2679        int iRow = permute_[i];
2680        regionX[iRow]=value;
2681      }
2682    }
2683    for (i=firstDense_-1;i>=0;i--) {
2684      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2685      CoinWorkDouble value=work[i]*diagonal_[i];
2686      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2687        int iRow = choleskyRow_[j+offset];
2688        value -= sparseFactor_[j]*work[iRow];
2689      }
2690      work[i]=value;
2691      int iRow = permute_[i];
2692      regionX[iRow]=value;
2693    }
2694    double largest=0.0;
2695    double largestV=0.0;
2696    for (i=0;i<numberRows_;i++) {
2697      largest = CoinMax(largest,CoinAbs(region[i]-regionX[i]));
2698      largestV = CoinMax(largestV,CoinAbs(region[i]));
2699    }
2700    printf("largest difference %g, largest %g, largest original %g\n",
2701           largest,largestV,largestO);
2702    delete [] regionX;
2703  }
2704#endif
2705}
2706#if CLP_LONG_CHOLESKY
2707/* Uses factorization to solve. */
2708void 
2709ClpCholeskyBase::solve (CoinWorkDouble * region) 
2710{
2711  assert (!whichDense_) ;
2712  CoinWorkDouble * work = reinterpret_cast<CoinWorkDouble *> (workDouble_);
2713  int i;
2714  CoinBigIndex j;
2715  for (i=0;i<numberRows_;i++) {
2716    int iRow = permute_[i];
2717    work[i] = region[iRow];
2718  }
2719  for (i=0;i<firstDense_;i++) {
2720    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2721    CoinWorkDouble value=work[i];
2722    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2723      int iRow = choleskyRow_[j+offset];
2724      work[iRow] -= sparseFactor_[j]*value;
2725    }
2726  }
2727  if (firstDense_<numberRows_) {
2728    // do dense
2729    ClpCholeskyDense dense;
2730    // just borrow space
2731    int nDense = numberRows_-firstDense_;
2732    dense.reserveSpace(this,nDense);
2733#if CLP_LONG_CHOLESKY!=1
2734    dense.solveLong(work+firstDense_);
2735#else
2736    dense.solveLongWork(work+firstDense_);
2737#endif
2738    for (i=numberRows_-1;i>=firstDense_;i--) {
2739      CoinWorkDouble value=work[i];
2740      int iRow = permute_[i];
2741      region[iRow]=value;
2742    }
2743  }
2744  for (i=firstDense_-1;i>=0;i--) {
2745    CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2746    CoinWorkDouble value=work[i]*diagonal_[i];
2747    for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2748      int iRow = choleskyRow_[j+offset];
2749      value -= sparseFactor_[j]*work[iRow];
2750    }
2751    work[i]=value;
2752    int iRow = permute_[i];
2753    region[iRow]=value;
2754  }
2755}
2756#endif
Note: See TracBrowser for help on using the repository browser.