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