source: branches/BSP/trunk/Clp/src/ClpCholeskyBase.cpp @ 1071

Last change on this file since 1071 was 894, checked in by forrest, 13 years ago

out messages

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 82.5 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_ = (double *) ClpCopyOfArray((longWork *) rhs.workDouble_,numberRows_);
89#endif
90  link_ = ClpCopyOfArray(rhs.link_,numberRows_);
91  workInteger_ = ClpCopyOfArray(rhs.workInteger_,numberRows_);
92  clique_ = ClpCopyOfArray(rhs.clique_,numberRows_);
93  memcpy(integerParameters_,rhs.integerParameters_,64*sizeof(int));
94  memcpy(doubleParameters_,rhs.doubleParameters_,64*sizeof(double));
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_ = (double *) ClpCopyOfArray((longWork *) 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 (double * region1, double * region2, const double * diagonal,
198                           double diagonalScaleFactor)
199{
200  if (!doKKT_) {
201    int iColumn;
202    int numberColumns = model_->numberColumns();
203    int numberTotal = numberRows_+numberColumns;
204    double * region1Save = new double[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    double maximumRHS = maximumAbsElement(region2,numberRows_);
212    double scale=1.0;
213    double unscale=1.0;
214    if (maximumRHS>1.0e-30) {
215      if (maximumRHS<=0.5) {
216        double 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        double 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    double * array = new double [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]&&fabs(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]&&fabs(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    memcpy(permute_+nn,used,np*sizeof(int));
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        memcpy(indices,choleskyRow_,sizeIndex_*sizeof(int));
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 += ((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_ = (double *) (new longWork[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 = (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 = (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 ((double) 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 double * 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  longDouble 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=sqrt(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  double largest=1.0;
1551  double smallest=COIN_DBL_MAX;
1552  int numberDense=0;
1553  if (!doKKT_) {
1554    const double * 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    longDouble delta2 = model_->delta(); // add delta*delta to diagonal
1580    delta2 *= delta2;
1581    // largest in initial matrix
1582    double 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            longDouble 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                longDouble value=element[j]*multiplier;
1605                work[jNewRow] += value;
1606              }
1607            }
1608          }
1609        }
1610        diagonal_[iRow]=work[iRow];
1611        largest2 = CoinMax(largest2,fabs(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,fabs(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,1.0e-11);
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        longDouble 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]=sqrt(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        solveLong(a,1);
1676      }
1677      dense_->resetRowsDropped();
1678      longDouble * denseBlob = dense_->aMatrix();
1679      longDouble * denseDiagonal = dense_->diagonal();
1680      // Update dense matrix
1681      for (i=0;i<numberDense;i++) {
1682        const longDouble * a = denseColumn_+i*numberRows_;
1683        // do diagonal
1684        longDouble value = denseDiagonal[i];
1685        const longDouble * b = denseColumn_+i*numberRows_;
1686        for (int k=0;k<numberRows_;k++) 
1687          value += a[k]*b[k];
1688        denseDiagonal[i]=value;
1689        for (int j=i+1;j<numberDense;j++) {
1690          longDouble value = 0.0;
1691          const longDouble * b = denseColumn_+j*numberRows_;
1692          for (int k=0;k<numberRows_;k++) 
1693            value += a[k]*b[k];
1694          *denseBlob=value;
1695          denseBlob++;
1696        }
1697      }
1698      // dense cholesky (? long double)
1699      int * dropped = new int [numberDense];
1700      dense_->factorizePart2(dropped);
1701      delete [] dropped;
1702    }
1703    // try allowing all every time
1704    //printf("trying ?\n");
1705    //for (iRow=0;iRow<numberRows_;iRow++) {
1706    //rowsDropped[iRow]=0;
1707    //rowsDropped_[iRow]=0;
1708    //}
1709    bool cleanCholesky;
1710    //if (model_->numberIterations()<20||(model_->numberIterations()&1)==0)
1711    if (model_->numberIterations()<2000) 
1712      cleanCholesky=true;
1713    else 
1714      cleanCholesky=false;
1715    if (cleanCholesky) {
1716      //drop fresh makes some formADAT easier
1717      if (newDropped||numberRowsDropped_) {
1718        newDropped=0;
1719        for (int i=0;i<numberRows_;i++) {
1720          char dropped = rowsDropped[i];
1721          rowsDropped_[i]=dropped;
1722          rowsDropped_[i]=0;
1723          if (dropped==2) {
1724            //dropped this time
1725            rowsDropped[newDropped++]=i;
1726            rowsDropped_[i]=0;
1727          } 
1728        } 
1729        numberRowsDropped_=newDropped;
1730        newDropped=-(2+newDropped);
1731      } 
1732    } else {
1733      if (newDropped) {
1734        newDropped=0;
1735        for (int i=0;i<numberRows_;i++) {
1736          char dropped = rowsDropped[i];
1737          rowsDropped_[i]=dropped;
1738          if (dropped==2) {
1739            //dropped this time
1740            rowsDropped[newDropped++]=i;
1741            rowsDropped_[i]=1;
1742          } 
1743        } 
1744      } 
1745      numberRowsDropped_+=newDropped;
1746      if (numberRowsDropped_&&0) {
1747        std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
1748          numberRowsDropped_<<" dropped)";
1749        if (newDropped) {
1750          std::cout<<" ( "<<newDropped<<" dropped this time)";
1751        } 
1752        std::cout<<std::endl;
1753      } 
1754    }
1755  } else {
1756    //KKT
1757    CoinPackedMatrix * quadratic = NULL;
1758    ClpQuadraticObjective * quadraticObj = 
1759      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
1760    if (quadraticObj) 
1761      quadratic = quadraticObj->quadraticObjective();
1762    // matrix
1763    int numberRowsModel = model_->numberRows();
1764    int numberColumns = model_->numberColumns();
1765    int numberTotal = numberColumns + numberRowsModel;
1766    // temp
1767    bool permuted=false;
1768    for (iRow=0;iRow<numberRows_;iRow++) {
1769      if (permute_[iRow]!=iRow) {
1770        permuted=true;
1771        break;
1772      }
1773    }
1774    // but fake it
1775    for (iRow=0;iRow<numberRows_;iRow++) {
1776      //permute_[iRow]=iRow; // force no permute
1777      //permuteInverse_[permute_[iRow]]=iRow;
1778    }
1779    if (permuted) {
1780      longDouble delta2 = model_->delta(); // add delta*delta to bottom
1781      delta2 *= delta2;
1782      // Need to permute - ugly
1783      if (!quadratic) {
1784        for (iRow=0;iRow<numberRows_;iRow++) {
1785          longDouble * put = sparseFactor_+choleskyStart_[iRow];
1786          int * which = choleskyRow_+indexStart_[iRow];
1787          int iOriginalRow = permute_[iRow];
1788          if (iOriginalRow<numberColumns) {
1789            iColumn=iOriginalRow;
1790            longDouble value = diagonal[iColumn];
1791            if (fabs(value)>1.0e-100) {
1792              value = 1.0/value;
1793              largest = CoinMax(largest,fabs(value));
1794              diagonal_[iRow] = -value;
1795              CoinBigIndex start=columnStart[iColumn];
1796              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1797              for (CoinBigIndex j=start;j<end;j++) {
1798                int kRow = row[j]+numberTotal;
1799                kRow=permuteInverse_[kRow];
1800                if (kRow>iRow) {
1801                  work[kRow]=element[j];
1802                  largest = CoinMax(largest,fabs(element[j]));
1803                }
1804              }
1805            } else {
1806              diagonal_[iRow] = -value;
1807            }
1808          } else if (iOriginalRow<numberTotal) {
1809            longDouble value = diagonal[iOriginalRow];
1810            if (fabs(value)>1.0e-100) {
1811              value = 1.0/value;
1812              largest = CoinMax(largest,fabs(value));
1813            } else {
1814              value = 1.0e100;
1815            }
1816            diagonal_[iRow] = -value;
1817            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1818            if (kRow>iRow) 
1819              work[kRow]=-1.0;
1820          } else {
1821            diagonal_[iRow]=delta2;
1822            int kRow = iOriginalRow-numberTotal;
1823            CoinBigIndex start=rowStart[kRow];
1824            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1825            for (CoinBigIndex j=start;j<end;j++) {
1826              int jRow = column[j];
1827              int jNewRow = permuteInverse_[jRow];
1828              if (jNewRow>iRow) {
1829                work[jNewRow]=elementByRow[j];
1830                largest = CoinMax(largest,fabs(elementByRow[j]));
1831              }
1832            }
1833            // slack - should it be permute
1834            kRow = permuteInverse_[kRow+numberColumns];
1835            if (kRow>iRow)
1836              work[kRow]=-1.0;
1837          }
1838          CoinBigIndex j;
1839          int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1840          for (j=0;j<number;j++) {
1841            int jRow = which[j];
1842            put[j]=work[jRow];
1843            work[jRow]=0.0;
1844          }
1845        }
1846      } else {
1847        // quadratic
1848        const int * columnQuadratic = quadratic->getIndices();
1849        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1850        const int * columnQuadraticLength = quadratic->getVectorLengths();
1851        const double * quadraticElement = quadratic->getElements();
1852        for (iRow=0;iRow<numberRows_;iRow++) {
1853          longDouble * put = sparseFactor_+choleskyStart_[iRow];
1854          int * which = choleskyRow_+indexStart_[iRow];
1855          int iOriginalRow = permute_[iRow];
1856          if (iOriginalRow<numberColumns) {
1857            CoinBigIndex j;
1858            iColumn=iOriginalRow;
1859            longDouble value = diagonal[iColumn];
1860            if (fabs(value)>1.0e-100) {
1861              value = 1.0/value;
1862              for (j=columnQuadraticStart[iColumn];
1863                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1864                int jColumn = columnQuadratic[j];
1865                int jNewColumn = permuteInverse_[jColumn];
1866                if (jNewColumn>iRow) {
1867                  work[jNewColumn]=-quadraticElement[j];
1868                } else if (iColumn==jColumn) {
1869                  value += quadraticElement[j];
1870                }
1871              }
1872              largest = CoinMax(largest,fabs(value));
1873              diagonal_[iRow] = -value;
1874              CoinBigIndex start=columnStart[iColumn];
1875              CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1876              for (j=start;j<end;j++) {
1877                int kRow = row[j]+numberTotal;
1878                kRow=permuteInverse_[kRow];
1879                if (kRow>iRow) {
1880                  work[kRow]=element[j];
1881                  largest = CoinMax(largest,fabs(element[j]));
1882                }
1883              }
1884            } else {
1885              diagonal_[iRow] = -value;
1886            }
1887          } else if (iOriginalRow<numberTotal) {
1888            longDouble value = diagonal[iOriginalRow];
1889            if (fabs(value)>1.0e-100) {
1890              value = 1.0/value;
1891              largest = CoinMax(largest,fabs(value));
1892            } else {
1893              value = 1.0e100;
1894            }
1895            diagonal_[iRow] = -value;
1896            int kRow = permuteInverse_[iOriginalRow+numberRowsModel];
1897            if (kRow>iRow) 
1898              work[kRow]=-1.0;
1899          } else {
1900            diagonal_[iRow]=delta2;
1901            int kRow = iOriginalRow-numberTotal;
1902            CoinBigIndex start=rowStart[kRow];
1903            CoinBigIndex end=rowStart[kRow]+rowLength[kRow];
1904            for (CoinBigIndex j=start;j<end;j++) {
1905              int jRow = column[j];
1906              int jNewRow = permuteInverse_[jRow];
1907              if (jNewRow>iRow) {
1908                work[jNewRow]=elementByRow[j];
1909                largest = CoinMax(largest,fabs(elementByRow[j]));
1910              }
1911            }
1912            // slack - should it be permute
1913            kRow = permuteInverse_[kRow+numberColumns];
1914            if (kRow>iRow)
1915              work[kRow]=-1.0;
1916          }
1917          CoinBigIndex j;
1918          int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
1919          for (j=0;j<number;j++) {
1920            int jRow = which[j];
1921            put[j]=work[jRow];
1922            work[jRow]=0.0;
1923          }
1924          for (j=0;j<numberRows_;j++)
1925            assert (!work[j]);
1926        }
1927      }
1928    } else {
1929      if (!quadratic) {
1930        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1931          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
1932          int * which = choleskyRow_+indexStart_[iColumn];
1933          longDouble value = diagonal[iColumn];
1934          if (fabs(value)>1.0e-100) {
1935            value = 1.0/value;
1936            largest = CoinMax(largest,fabs(value));
1937            diagonal_[iColumn] = -value;
1938            CoinBigIndex start=columnStart[iColumn];
1939            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1940            for (CoinBigIndex j=start;j<end;j++) {
1941              //choleskyRow_[numberElements]=row[j]+numberTotal;
1942              //sparseFactor_[numberElements++]=element[j];
1943              work[row[j]+numberTotal]=element[j];
1944              largest = CoinMax(largest,fabs(element[j]));
1945            }
1946          } else {
1947            diagonal_[iColumn] = -value;
1948          }
1949          CoinBigIndex j;
1950          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
1951          for (j=0;j<number;j++) {
1952            int jRow = which[j];
1953            put[j]=work[jRow];
1954            work[jRow]=0.0;
1955          }
1956        }
1957      } else {
1958        // Quadratic
1959        const int * columnQuadratic = quadratic->getIndices();
1960        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1961        const int * columnQuadraticLength = quadratic->getVectorLengths();
1962        const double * quadraticElement = quadratic->getElements();
1963        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1964          longDouble * put = sparseFactor_+choleskyStart_[iColumn];
1965          int * which = choleskyRow_+indexStart_[iColumn];
1966          int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
1967          longDouble value = diagonal[iColumn];
1968          CoinBigIndex j;
1969          if (fabs(value)>1.0e-100) {
1970            value = 1.0/value;
1971            for (j=columnQuadraticStart[iColumn];
1972                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1973              int jColumn = columnQuadratic[j];
1974              if (jColumn>iColumn) {
1975                work[jColumn]=-quadraticElement[j];
1976              } else if (iColumn==jColumn) {
1977                value += quadraticElement[j];
1978              }
1979            }
1980            largest = CoinMax(largest,fabs(value));
1981            diagonal_[iColumn] = -value;
1982            CoinBigIndex start=columnStart[iColumn];
1983            CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
1984            for (j=start;j<end;j++) {
1985              work[row[j]+numberTotal]=element[j];
1986              largest = CoinMax(largest,fabs(element[j]));
1987            }
1988            for (j=0;j<number;j++) {
1989              int jRow = which[j];
1990              put[j]=work[jRow];
1991              work[jRow]=0.0;
1992            }
1993          } else {
1994            value = 1.0e100;
1995            diagonal_[iColumn] = -value;
1996            for (j=0;j<number;j++) {
1997              int jRow = which[j];
1998              put[j]=work[jRow];
1999            }
2000          }
2001        }
2002      }
2003      // slacks
2004      for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
2005        longDouble * put = sparseFactor_+choleskyStart_[iColumn];
2006        int * which = choleskyRow_+indexStart_[iColumn];
2007        longDouble value = diagonal[iColumn];
2008        if (fabs(value)>1.0e-100) {
2009          value = 1.0/value;
2010          largest = CoinMax(largest,fabs(value));
2011        } else {
2012          value = 1.0e100;
2013        }
2014        diagonal_[iColumn] = -value;
2015        work[iColumn-numberColumns+numberTotal]=-1.0;
2016        CoinBigIndex j;
2017        int number = choleskyStart_[iColumn+1]-choleskyStart_[iColumn];
2018        for (j=0;j<number;j++) {
2019          int jRow = which[j];
2020          put[j]=work[jRow];
2021          work[jRow]=0.0;
2022        }
2023      }
2024      // Finish diagonal
2025      longDouble delta2 = model_->delta(); // add delta*delta to bottom
2026      delta2 *= delta2;
2027      for (iRow=0;iRow<numberRowsModel;iRow++) {
2028        longDouble * put = sparseFactor_+choleskyStart_[iRow+numberTotal];
2029        diagonal_[iRow+numberTotal]=delta2;
2030        CoinBigIndex j;
2031        int number = choleskyStart_[iRow+numberTotal+1]-choleskyStart_[iRow+numberTotal];
2032        for (j=0;j<number;j++) {
2033          put[j]=0.0;
2034        }
2035      }
2036    }
2037    //check sizes
2038    largest*=1.0e-20;
2039    largest = CoinMin(largest,1.0e-11);
2040    doubleParameters_[10]=CoinMax(1.0e-20,largest);
2041    integerParameters_[20]=0;
2042    doubleParameters_[3]=0.0;
2043    doubleParameters_[4]=COIN_DBL_MAX;
2044    // Set up LDL cutoff
2045    integerParameters_[34]=numberTotal;
2046    // KKT
2047    int * rowsDropped2 = new int[numberRows_];
2048    CoinZeroN(rowsDropped2,numberRows_);
2049    factorizePart2(rowsDropped2);
2050    newDropped=integerParameters_[20];
2051    largest=doubleParameters_[3];
2052    smallest=doubleParameters_[4];
2053    if (model_->messageHandler()->logLevel()>1) 
2054      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
2055    choleskyCondition_=largest/smallest;
2056    // Should save adjustments in ..R_
2057    int n1=0,n2=0;
2058    double * primalR = model_->primalR();
2059    double * dualR = model_->dualR();
2060    for (iRow=0;iRow<numberTotal;iRow++) { 
2061      if (rowsDropped2[iRow]) {
2062        n1++;
2063        //printf("row region1 %d dropped\n",iRow);
2064        //rowsDropped_[iRow]=1;
2065        rowsDropped_[iRow]=0;
2066        primalR[iRow]=doubleParameters_[20];
2067      } else {
2068        rowsDropped_[iRow]=0;
2069        primalR[iRow]=0.0;
2070      }
2071    }
2072    for (;iRow<numberRows_;iRow++) {
2073      if (rowsDropped2[iRow]) {
2074        n2++;
2075        //printf("row region2 %d dropped\n",iRow);
2076        //rowsDropped_[iRow]=1;
2077        rowsDropped_[iRow]=0;
2078        dualR[iRow-numberTotal]=doubleParameters_[34];
2079      } else {
2080        rowsDropped_[iRow]=0;
2081        dualR[iRow-numberTotal]=0.0;
2082      }
2083    }
2084    delete [] rowsDropped2;
2085  }
2086  status_=0;
2087  return newDropped;
2088}
2089/* Factorize - filling in rowsDropped and returning number dropped
2090   in integerParam.
2091*/
2092void 
2093ClpCholeskyBase::factorizePart2(int * rowsDropped) 
2094{
2095  longDouble largest=doubleParameters_[3];
2096  longDouble smallest=doubleParameters_[4];
2097  int numberDropped=integerParameters_[20];
2098  // probably done before
2099  largest=0.0;
2100  smallest=COIN_DBL_MAX;
2101  numberDropped=0;
2102  double dropValue = doubleParameters_[10];
2103  int firstPositive=integerParameters_[34];
2104  longDouble * d = ClpCopyOfArray(diagonal_,numberRows_);
2105  int iRow;
2106  // minimum size before clique done
2107  //#define MINCLIQUE INT_MAX
2108#define MINCLIQUE 3
2109  longDouble * work = workDouble_;
2110  CoinBigIndex * first = workInteger_;
2111 
2112  for (iRow=0;iRow<numberRows_;iRow++) {
2113    link_[iRow]=-1;
2114    work[iRow]=0.0;
2115    first[iRow]=choleskyStart_[iRow];
2116  }
2117
2118  int lastClique=-1;
2119  bool inClique=false;
2120  bool newClique=false;
2121  bool endClique=false;
2122  int lastRow=0;
2123  int cliqueSize=0;
2124  CoinBigIndex cliquePointer=0;
2125  int nextRow2=-1;
2126 
2127  for (iRow=0;iRow<firstDense_+1;iRow++) {
2128    if (iRow<firstDense_) {
2129      endClique=false;
2130      if (clique_[iRow]>0) {
2131        // this is in a clique
2132        inClique=true;
2133        if (clique_[iRow]>lastClique) {
2134          // new Clique
2135          newClique=true;
2136          // If we have clique going then signal to do old one
2137          endClique=(lastClique>0);
2138        } else {
2139          // Still in clique
2140          newClique=false;
2141        }
2142      } else {
2143        // not in clique
2144        inClique=false;
2145        newClique=false;
2146        // If we have clique going then signal to do old one
2147        endClique=(lastClique>0);
2148      }
2149      lastClique=clique_[iRow];
2150    } else if (inClique) {
2151      // Finish off
2152      endClique=true;
2153    } else {
2154      break;
2155    }
2156    if (endClique) {
2157      // We have just finished updating a clique - do block pivot and clean up
2158      int jRow;
2159      for ( jRow=lastRow;jRow<iRow;jRow++) {
2160        int jCount = jRow-lastRow;
2161        longDouble diagonalValue = diagonal_[jRow];
2162        CoinBigIndex start=choleskyStart_[jRow];
2163        CoinBigIndex end=choleskyStart_[jRow+1];
2164        for (int kRow=lastRow;kRow<jRow;kRow++) {
2165          jCount--;
2166          CoinBigIndex get = choleskyStart_[kRow]+jCount;
2167          longDouble a_jk = sparseFactor_[get];
2168          longDouble value1 = d[kRow]*a_jk;
2169          diagonalValue -= a_jk*value1;
2170          for (CoinBigIndex j=start;j<end;j++)
2171            sparseFactor_[j] -= value1*sparseFactor_[++get];
2172        }
2173        // check
2174        int originalRow = permute_[jRow];
2175        if (originalRow<firstPositive) {
2176          // must be negative
2177          if (diagonalValue<=-dropValue) {
2178            smallest = CoinMin(smallest,-diagonalValue);
2179            largest = CoinMax(largest,-diagonalValue);
2180            d[jRow]=diagonalValue;
2181            diagonalValue = 1.0/diagonalValue;
2182          } else {
2183            rowsDropped[originalRow]=2;
2184            d[jRow]=-1.0e100;
2185            diagonalValue=0.0;
2186            integerParameters_[20]++;
2187          }
2188        } else {
2189          // must be positive
2190          if (diagonalValue>=dropValue) {
2191            smallest = CoinMin(smallest,diagonalValue);
2192            largest = CoinMax(largest,diagonalValue);
2193            d[jRow]=diagonalValue;
2194            diagonalValue = 1.0/diagonalValue;
2195          } else {
2196            rowsDropped[originalRow]=2;
2197            d[jRow]=1.0e100;
2198            diagonalValue=0.0;
2199            integerParameters_[20]++;
2200          }
2201        }
2202        diagonal_[jRow]=diagonalValue;
2203        for (CoinBigIndex j=start;j<end;j++) {
2204          sparseFactor_[j] *= diagonalValue;
2205        }
2206      }
2207      if (nextRow2>=0) {
2208        for ( jRow=lastRow;jRow<iRow-1;jRow++) {
2209          link_[jRow]=jRow+1;
2210        }
2211        link_[iRow-1]=link_[nextRow2];
2212        link_[nextRow2]=lastRow;
2213      }
2214    }
2215    if (iRow==firstDense_)
2216      break; // we were just cleaning up
2217    if (newClique) {
2218      // initialize new clique
2219      lastRow=iRow;
2220      cliquePointer=choleskyStart_[iRow];
2221      cliqueSize=choleskyStart_[iRow+1]-cliquePointer+1;
2222    }
2223    // for each column L[*,kRow] that affects L[*,iRow]
2224    longDouble diagonalValue=diagonal_[iRow];
2225    int nextRow = link_[iRow];
2226    int kRow=0;
2227    while (1) {
2228      kRow=nextRow;
2229      if (kRow<0)
2230        break; // out of loop
2231      nextRow=link_[kRow];
2232      // Modify by outer product of L[*,irow] by L[*,krow] from first
2233      CoinBigIndex k=first[kRow];
2234      CoinBigIndex end=choleskyStart_[kRow+1];
2235      assert(k<end);
2236      longDouble a_ik=sparseFactor_[k++];
2237      longDouble value1=d[kRow]*a_ik;
2238      // update first
2239      first[kRow]=k;
2240      diagonalValue -= value1*a_ik;
2241      CoinBigIndex offset = indexStart_[kRow]-choleskyStart_[kRow];
2242      if (k<end) {
2243        int jRow = choleskyRow_[k+offset];
2244        if (clique_[kRow]<MINCLIQUE) {
2245          link_[kRow]=link_[jRow];
2246          link_[jRow]=kRow;
2247          for (;k<end;k++) {
2248            int jRow = choleskyRow_[k+offset];
2249            work[jRow] += sparseFactor_[k]*value1;
2250          }
2251        } else {
2252          // Clique
2253          CoinBigIndex currentIndex = k+offset;
2254          int linkSave=link_[jRow];
2255          link_[jRow]=kRow;
2256          work[kRow]=value1; // ? or a_jk
2257          int last = kRow+clique_[kRow];
2258          for (int kkRow=kRow+1;kkRow<last;kkRow++) {
2259            CoinBigIndex j=first[kkRow];
2260            //int iiRow = choleskyRow_[j+indexStart_[kkRow]-choleskyStart_[kkRow]];
2261            longDouble a = sparseFactor_[j];
2262            longDouble dValue = d[kkRow]*a;
2263            diagonalValue -= a*dValue;
2264            work[kkRow]=dValue;
2265            first[kkRow]++;
2266            link_[kkRow-1]=kkRow;
2267          }
2268          nextRow = link_[last-1];
2269          link_[last-1]=linkSave;
2270          int length = end-k;
2271          for (int i=0;i<length;i++) {
2272            int lRow = choleskyRow_[currentIndex++];
2273            longDouble t0 = work[lRow];
2274            for (int kkRow=kRow;kkRow<last;kkRow++) {
2275              CoinBigIndex j = first[kkRow]+i;
2276              t0 += work[kkRow]*sparseFactor_[j];
2277            }
2278            work[lRow]=t0;
2279          }
2280        }
2281      }
2282    }
2283    // Now apply
2284    if (inClique) {
2285      // in clique
2286      diagonal_[iRow]=diagonalValue;
2287      CoinBigIndex start=choleskyStart_[iRow];
2288      CoinBigIndex end=choleskyStart_[iRow+1];
2289      CoinBigIndex currentIndex = indexStart_[iRow];
2290      nextRow2=-1;
2291      CoinBigIndex get=start+clique_[iRow]-1;
2292      if (get<end) {
2293        nextRow2 = choleskyRow_[currentIndex+get-start];
2294        first[iRow]=get;
2295      }
2296      for (CoinBigIndex j=start;j<end;j++) {
2297        int kRow = choleskyRow_[currentIndex++];
2298        sparseFactor_[j] -= work[kRow]; // times?
2299        work[kRow]=0.0;
2300      }
2301    } else {
2302      // not in clique
2303      int originalRow = permute_[iRow];
2304      if (originalRow<firstPositive) {
2305        // must be negative
2306        if (diagonalValue<=-dropValue) {
2307          smallest = CoinMin(smallest,-diagonalValue);
2308          largest = CoinMax(largest,-diagonalValue);
2309          d[iRow]=diagonalValue;
2310          diagonalValue = 1.0/diagonalValue;
2311        } else {
2312          rowsDropped[originalRow]=2;
2313          d[iRow]=-1.0e100;
2314          diagonalValue=0.0;
2315          integerParameters_[20]++;
2316        }
2317      } else {
2318        // must be positive
2319        if (diagonalValue>=dropValue) {
2320          smallest = CoinMin(smallest,diagonalValue);
2321          largest = CoinMax(largest,diagonalValue);
2322          d[iRow]=diagonalValue;
2323          diagonalValue = 1.0/diagonalValue;
2324        } else {
2325          rowsDropped[originalRow]=2;
2326          d[iRow]=1.0e100;
2327          diagonalValue=0.0;
2328          integerParameters_[20]++;
2329        }
2330      }
2331      diagonal_[iRow]=diagonalValue;
2332      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
2333      CoinBigIndex start = choleskyStart_[iRow];
2334      CoinBigIndex end = choleskyStart_[iRow+1];
2335      assert (first[iRow]==start);
2336      if (start<end) {
2337        int nextRow = choleskyRow_[start+offset];
2338        link_[iRow]=link_[nextRow];
2339        link_[nextRow]=iRow;
2340        for (CoinBigIndex j=start;j<end;j++) {
2341          int jRow = choleskyRow_[j+offset];
2342          longDouble value = sparseFactor_[j]-work[jRow];
2343          work[jRow]=0.0;
2344          sparseFactor_[j]= diagonalValue*value;
2345        }
2346      }
2347    }
2348  }
2349  if (firstDense_<numberRows_) {
2350    // do dense
2351    // update dense part
2352    updateDense(d,work,first);
2353    ClpCholeskyDense dense;
2354    // just borrow space
2355    int nDense = numberRows_-firstDense_;
2356    if (doKKT_) {
2357      for (iRow=firstDense_;iRow<numberRows_;iRow++) {
2358        int originalRow=permute_[iRow];
2359        if (originalRow>=firstPositive) {
2360          firstPositive=iRow-firstDense_;
2361          break;
2362        }
2363      }
2364    }
2365    dense.reserveSpace(this,nDense);
2366    int * dropped = new int[nDense];
2367    memset(dropped,0,nDense*sizeof(int));
2368    dense.setDoubleParameter(3,largest);
2369    dense.setDoubleParameter(4,smallest);
2370    dense.setDoubleParameter(10,dropValue);
2371    dense.setIntegerParameter(20,0);
2372    dense.setIntegerParameter(34,firstPositive);
2373    dense.factorizePart2(dropped);
2374    largest=dense.getDoubleParameter(3);
2375    smallest=dense.getDoubleParameter(4);
2376    integerParameters_[20]+=dense.getIntegerParameter(20);
2377    for (iRow=firstDense_;iRow<numberRows_;iRow++) {
2378      int originalRow=permute_[iRow];
2379      rowsDropped[originalRow]=dropped[iRow-firstDense_];
2380    }
2381    delete [] dropped;
2382  }
2383  delete [] d;
2384  doubleParameters_[3]=largest;
2385  doubleParameters_[4]=smallest;
2386  return;
2387}
2388// Updates dense part (broken out for profiling)
2389void ClpCholeskyBase::updateDense(longDouble * d, longDouble * work, int * first)
2390{
2391  for (int iRow=0;iRow<firstDense_;iRow++) {
2392    CoinBigIndex start=first[iRow];
2393    CoinBigIndex end=choleskyStart_[iRow+1];
2394    if (start<end) {
2395      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
2396      if (clique_[iRow]<2) {
2397        longDouble dValue=d[iRow];
2398        for (CoinBigIndex k=start;k<end;k++) {
2399          int kRow = choleskyRow_[k+offset]; 
2400          assert(kRow>=firstDense_);
2401          longDouble a_ik=sparseFactor_[k];
2402          longDouble value1=dValue*a_ik;
2403          diagonal_[kRow] -= value1*a_ik;
2404          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2405          for (CoinBigIndex j=k+1;j<end;j++) {
2406            int jRow=choleskyRow_[j+offset];
2407            longDouble a_jk = sparseFactor_[j];
2408            sparseFactor_[base+jRow] -= a_jk*value1;
2409          }
2410        }
2411      } else if (clique_[iRow]<3) {
2412        // do as pair
2413        longDouble dValue0=d[iRow];
2414        longDouble dValue1=d[iRow+1];
2415        int offset1 = first[iRow+1]-start;
2416        // skip row
2417        iRow++;
2418        for (CoinBigIndex k=start;k<end;k++) {
2419          int kRow = choleskyRow_[k+offset]; 
2420          assert(kRow>=firstDense_);
2421          longDouble a_ik0=sparseFactor_[k];
2422          longDouble value0=dValue0*a_ik0;
2423          longDouble a_ik1=sparseFactor_[k+offset1];
2424          longDouble value1=dValue1*a_ik1;
2425          diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1;
2426          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2427          for (CoinBigIndex j=k+1;j<end;j++) {
2428            int jRow=choleskyRow_[j+offset];
2429            longDouble a_jk0 = sparseFactor_[j];
2430            longDouble a_jk1 = sparseFactor_[j+offset1];
2431            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1;
2432          }
2433        }
2434#define MANY_REGISTERS
2435#ifdef MANY_REGISTERS
2436      } else if (clique_[iRow]==3) {
2437#else
2438      } else {
2439#endif
2440        // do as clique
2441        // maybe later get fancy on big cliques and do transpose copy
2442        // seems only worth going to 3 on Intel
2443        longDouble dValue0=d[iRow];
2444        longDouble dValue1=d[iRow+1];
2445        longDouble dValue2=d[iRow+2];
2446        // get offsets and skip rows
2447        int offset1 = first[++iRow]-start;
2448        int offset2 = first[++iRow]-start;
2449        for (CoinBigIndex k=start;k<end;k++) {
2450          int kRow = choleskyRow_[k+offset]; 
2451          assert(kRow>=firstDense_);
2452          double diagonalValue=diagonal_[kRow];
2453          longDouble a_ik0=sparseFactor_[k];
2454          longDouble value0=dValue0*a_ik0;
2455          longDouble a_ik1=sparseFactor_[k+offset1];
2456          longDouble value1=dValue1*a_ik1;
2457          longDouble a_ik2=sparseFactor_[k+offset2];
2458          longDouble value2=dValue2*a_ik2;
2459          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2460          diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2;
2461          for (CoinBigIndex j=k+1;j<end;j++) {
2462            int jRow=choleskyRow_[j+offset];
2463            longDouble a_jk0 = sparseFactor_[j];
2464            longDouble a_jk1 = sparseFactor_[j+offset1];
2465            longDouble a_jk2 = sparseFactor_[j+offset2];
2466            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2;
2467          }
2468        }
2469#ifdef MANY_REGISTERS
2470      } else {
2471        // do as clique
2472        // maybe later get fancy on big cliques and do transpose copy
2473        // maybe only worth going to 3 on Intel (but may have hidden registers)
2474        longDouble dValue0=d[iRow];
2475        longDouble dValue1=d[iRow+1];
2476        longDouble dValue2=d[iRow+2];
2477        longDouble dValue3=d[iRow+3];
2478        // get offsets and skip rows
2479        int offset1 = first[++iRow]-start;
2480        int offset2 = first[++iRow]-start;
2481        int offset3 = first[++iRow]-start;
2482        for (CoinBigIndex k=start;k<end;k++) {
2483          int kRow = choleskyRow_[k+offset]; 
2484          assert(kRow>=firstDense_);
2485          double diagonalValue=diagonal_[kRow];
2486          longDouble a_ik0=sparseFactor_[k];
2487          longDouble value0=dValue0*a_ik0;
2488          longDouble a_ik1=sparseFactor_[k+offset1];
2489          longDouble value1=dValue1*a_ik1;
2490          longDouble a_ik2=sparseFactor_[k+offset2];
2491          longDouble value2=dValue2*a_ik2;
2492          longDouble a_ik3=sparseFactor_[k+offset3];
2493          longDouble value3=dValue3*a_ik3;
2494          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
2495          diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3);
2496          for (CoinBigIndex j=k+1;j<end;j++) {
2497            int jRow=choleskyRow_[j+offset];
2498            longDouble a_jk0 = sparseFactor_[j];
2499            longDouble a_jk1 = sparseFactor_[j+offset1];
2500            longDouble a_jk2 = sparseFactor_[j+offset2];
2501            longDouble a_jk3 = sparseFactor_[j+offset3];
2502            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3;
2503          }
2504        }
2505#endif
2506      }
2507    }
2508  }
2509}
2510/* Uses factorization to solve. */
2511void 
2512ClpCholeskyBase::solve (double * region) 
2513{
2514  if (!whichDense_) {
2515    solve(region,3);
2516  } else {
2517    // dense columns
2518    int i;
2519    solve(region,1);
2520    // do change;
2521    int numberDense = dense_->numberRows();
2522    double * change = new double[numberDense];
2523    for (i=0;i<numberDense;i++) {
2524      const longDouble * a = denseColumn_+i*numberRows_;
2525      longDouble value =0.0;
2526      for (int iRow=0;iRow<numberRows_;iRow++) 
2527        value += a[iRow]*region[iRow];
2528      change[i]=value;
2529    }
2530    // solve
2531    dense_->solve(change);
2532    for (i=0;i<numberDense;i++) {
2533      const longDouble * a = denseColumn_+i*numberRows_;
2534      longDouble value = change[i];
2535      for (int iRow=0;iRow<numberRows_;iRow++) 
2536        region[iRow] -= value*a[iRow];
2537    }
2538    delete [] change;
2539    // and finish off
2540    solve(region,2);
2541  }
2542}
2543/* solve - 1 just first half, 2 just second half - 3 both.
2544   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
2545*/
2546void 
2547ClpCholeskyBase::solve(double * region, int type)
2548{
2549#ifdef CLP_DEBUG
2550  double * regionX=NULL;
2551  if (sizeof(longWork)!=sizeof(double)&&type==3) {
2552    regionX=ClpCopyOfArray(region,numberRows_);
2553  }
2554#endif
2555  longWork * work = (longWork *) workDouble_;
2556  int i;
2557  CoinBigIndex j;
2558  for (i=0;i<numberRows_;i++) {
2559    int iRow = permute_[i];
2560    work[i] = region[iRow];
2561  }
2562  switch (type) {
2563  case 1:
2564    for (i=0;i<numberRows_;i++) {
2565      longDouble value=work[i];
2566      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2567      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2568        int iRow = choleskyRow_[j+offset];
2569        work[iRow] -= sparseFactor_[j]*value;
2570      }
2571    }
2572    for (i=0;i<numberRows_;i++) {
2573      int iRow = permute_[i];
2574      region[iRow]=work[i]*diagonal_[i];
2575    }
2576    break;
2577  case 2:
2578    for (i=numberRows_-1;i>=0;i--) {
2579      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2580      longDouble value=work[i]*diagonal_[i];
2581      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2582        int iRow = choleskyRow_[j+offset];
2583        value -= sparseFactor_[j]*work[iRow];
2584      }
2585      work[i]=value;
2586      int iRow = permute_[i];
2587      region[iRow]=value;
2588    }
2589    break;
2590  case 3:
2591    for (i=0;i<firstDense_;i++) {
2592      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2593      longDouble value=work[i];
2594      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2595        int iRow = choleskyRow_[j+offset];
2596        work[iRow] -= sparseFactor_[j]*value;
2597      }
2598    }
2599    if (firstDense_<numberRows_) {
2600      // do dense
2601      ClpCholeskyDense dense;
2602      // just borrow space
2603      int nDense = numberRows_-firstDense_;
2604      dense.reserveSpace(this,nDense);
2605#if CLP_LONG_CHOLESKY!=1
2606      dense.solveLong(work+firstDense_);
2607#else
2608      dense.solveLongWork(work+firstDense_);
2609#endif
2610      for (i=numberRows_-1;i>=firstDense_;i--) {
2611        longDouble value=work[i];
2612        int iRow = permute_[i];
2613        region[iRow]=value;
2614      }
2615    }
2616    for (i=firstDense_-1;i>=0;i--) {
2617      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2618      longDouble value=work[i]*diagonal_[i];
2619      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2620        int iRow = choleskyRow_[j+offset];
2621        value -= sparseFactor_[j]*work[iRow];
2622      }
2623      work[i]=value;
2624      int iRow = permute_[i];
2625      region[iRow]=value;
2626    }
2627    break;
2628  }
2629#ifdef CLP_DEBUG
2630  if (regionX) {
2631    longDouble * work = workDouble_;
2632    int i;
2633    CoinBigIndex j;
2634    double largestO=0.0;
2635    for (i=0;i<numberRows_;i++) {
2636      largestO = CoinMax(largestO,fabs(regionX[i]));
2637    }
2638    for (i=0;i<numberRows_;i++) {
2639      int iRow = permute_[i];
2640      work[i] = regionX[iRow];
2641    }
2642    for (i=0;i<firstDense_;i++) {
2643      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2644      longDouble value=work[i];
2645      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2646        int iRow = choleskyRow_[j+offset];
2647        work[iRow] -= sparseFactor_[j]*value;
2648      }
2649    }
2650    if (firstDense_<numberRows_) {
2651      // do dense
2652      ClpCholeskyDense dense;
2653      // just borrow space
2654      int nDense = numberRows_-firstDense_;
2655      dense.reserveSpace(this,nDense);
2656      dense.solveLong(work+firstDense_);
2657      for (i=numberRows_-1;i>=firstDense_;i--) {
2658        longDouble value=work[i];
2659        int iRow = permute_[i];
2660        regionX[iRow]=value;
2661      }
2662    }
2663    for (i=firstDense_-1;i>=0;i--) {
2664      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2665      longDouble value=work[i]*diagonal_[i];
2666      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2667        int iRow = choleskyRow_[j+offset];
2668        value -= sparseFactor_[j]*work[iRow];
2669      }
2670      work[i]=value;
2671      int iRow = permute_[i];
2672      regionX[iRow]=value;
2673    }
2674    double largest=0.0;
2675    double largestV=0.0;
2676    for (i=0;i<numberRows_;i++) {
2677      largest = CoinMax(largest,fabs(region[i]-regionX[i]));
2678      largestV = CoinMax(largestV,fabs(region[i]));
2679    }
2680    printf("largest difference %g, largest %g, largest original %g\n",
2681           largest,largestV,largestO);
2682    delete [] regionX;
2683  }
2684#endif
2685}
2686/* solve - 1 just first half, 2 just second half - 3 both.
2687   If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
2688*/
2689void 
2690ClpCholeskyBase::solveLong(longDouble * region, int type)
2691{
2692  int i;
2693  CoinBigIndex j;
2694  for (i=0;i<numberRows_;i++) {
2695    int iRow = permute_[i];
2696    workDouble_[i] = region[iRow];
2697  }
2698  switch (type) {
2699  case 1:
2700    for (i=0;i<numberRows_;i++) {
2701      longDouble value=workDouble_[i];
2702      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2703      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2704        int iRow = choleskyRow_[j+offset];
2705        workDouble_[iRow] -= sparseFactor_[j]*value;
2706      }
2707    }
2708    for (i=0;i<numberRows_;i++) {
2709      int iRow = permute_[i];
2710      region[iRow]=workDouble_[i]*diagonal_[i];
2711    }
2712    break;
2713  case 2:
2714    for (i=numberRows_-1;i>=0;i--) {
2715      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2716      longDouble value=workDouble_[i]*diagonal_[i];
2717      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2718        int iRow = choleskyRow_[j+offset];
2719        value -= sparseFactor_[j]*workDouble_[iRow];
2720      }
2721      workDouble_[i]=value;
2722      int iRow = permute_[i];
2723      region[iRow]=value;
2724    }
2725    break;
2726  case 3:
2727    for (i=0;i<firstDense_;i++) {
2728      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2729      longDouble value=workDouble_[i];
2730      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2731        int iRow = choleskyRow_[j+offset];
2732        workDouble_[iRow] -= sparseFactor_[j]*value;
2733      }
2734    }
2735    if (firstDense_<numberRows_) {
2736      // do dense
2737      ClpCholeskyDense dense;
2738      // just borrow space
2739      int nDense = numberRows_-firstDense_;
2740      dense.reserveSpace(this,nDense);
2741      dense.solveLong(workDouble_+firstDense_);
2742      for (i=numberRows_-1;i>=firstDense_;i--) {
2743        longDouble value=workDouble_[i];
2744        int iRow = permute_[i];
2745        region[iRow]=value;
2746      }
2747    }
2748    for (i=firstDense_-1;i>=0;i--) {
2749      CoinBigIndex offset = indexStart_[i]-choleskyStart_[i];
2750      longDouble value=workDouble_[i]*diagonal_[i];
2751      for (j=choleskyStart_[i];j<choleskyStart_[i+1];j++) {
2752        int iRow = choleskyRow_[j+offset];
2753        value -= sparseFactor_[j]*workDouble_[iRow];
2754      }
2755      workDouble_[i]=value;
2756      int iRow = permute_[i];
2757      region[iRow]=value;
2758    }
2759    break;
2760  }
2761}
2762
Note: See TracBrowser for help on using the repository browser.