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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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