source: trunk/Clp/src/ClpCholeskyUfl.cpp @ 2449

Last change on this file since 2449 was 2385, checked in by unxusr, 11 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 13.8 KB
Line 
1/* $Id: ClpCholeskyUfl.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2004, 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#include "ClpConfig.h"
7
8extern "C" {
9#ifndef COIN_HAS_CHOLMOD
10#ifndef COIN_HAS_AMD
11#error "Need to have AMD or CHOLMOD to compile ClpCholeskyUfl."
12#else
13#include "amd.h"
14#endif
15#else
16#include "cholmod.h"
17#endif
18}
19
20#include "CoinPragma.hpp"
21#include "ClpCholeskyUfl.hpp"
22#include "ClpMessage.hpp"
23#include "ClpInterior.hpp"
24#include "CoinHelperFunctions.hpp"
25#include "ClpHelperFunctions.hpp"
26//#############################################################################
27// Constructors / Destructor / Assignment
28//#############################################################################
29
30//-------------------------------------------------------------------
31// Default Constructor
32//-------------------------------------------------------------------
33ClpCholeskyUfl::ClpCholeskyUfl(int denseThreshold)
34  : ClpCholeskyBase(denseThreshold)
35{
36  type_ = 14;
37  L_ = NULL;
38  c_ = NULL;
39
40#ifdef COIN_HAS_CHOLMOD
41  c_ = (cholmod_common *)malloc(sizeof(cholmod_common));
42  cholmod_start(c_);
43  // Can't use supernodal as may not be positive definite
44  c_->supernodal = 0;
45#endif
46}
47
48//-------------------------------------------------------------------
49// Copy constructor
50//-------------------------------------------------------------------
51ClpCholeskyUfl::ClpCholeskyUfl(const ClpCholeskyUfl &rhs)
52  : ClpCholeskyBase(rhs)
53{
54  abort();
55}
56
57//-------------------------------------------------------------------
58// Destructor
59//-------------------------------------------------------------------
60ClpCholeskyUfl::~ClpCholeskyUfl()
61{
62#ifdef COIN_HAS_CHOLMOD
63  cholmod_free_factor(&L_, c_);
64  cholmod_finish(c_);
65  free(c_);
66#endif
67}
68
69//----------------------------------------------------------------
70// Assignment operator
71//-------------------------------------------------------------------
72ClpCholeskyUfl &
73ClpCholeskyUfl::operator=(const ClpCholeskyUfl &rhs)
74{
75  if (this != &rhs) {
76    ClpCholeskyBase::operator=(rhs);
77    abort();
78  }
79  return *this;
80}
81
82//-------------------------------------------------------------------
83// Clone
84//-------------------------------------------------------------------
85ClpCholeskyBase *ClpCholeskyUfl::clone() const
86{
87  return new ClpCholeskyUfl(*this);
88}
89
90#ifndef COIN_HAS_CHOLMOD
91/* Orders rows and saves pointer to matrix.and model */
92int ClpCholeskyUfl::order(ClpInterior *model)
93{
94  int iRow;
95  model_ = model;
96  if (preOrder(false, true, doKKT_))
97    return -1;
98  permuteInverse_ = new int[numberRows_];
99  permute_ = new int[numberRows_];
100  double Control[AMD_CONTROL];
101  double Info[AMD_INFO];
102
103  amd_defaults(Control);
104  //amd_control(Control);
105
106  int returnCode = amd_order(numberRows_, choleskyStart_, choleskyRow_,
107    permute_, Control, Info);
108  delete[] choleskyRow_;
109  choleskyRow_ = NULL;
110  delete[] choleskyStart_;
111  choleskyStart_ = NULL;
112  //amd_info(Info);
113
114  if (returnCode != AMD_OK) {
115    std::cout << "AMD ordering failed" << std::endl;
116    return 1;
117  }
118  for (iRow = 0; iRow < numberRows_; iRow++) {
119    permuteInverse_[permute_[iRow]] = iRow;
120  }
121  return 0;
122}
123#else
124/* Orders rows and saves pointer to matrix.and model */
125int ClpCholeskyUfl::order(ClpInterior *model)
126{
127  numberRows_ = model->numberRows();
128  if (doKKT_) {
129    numberRows_ += numberRows_ + model->numberColumns();
130    printf("finish coding UFL KKT!\n");
131    abort();
132  }
133  rowsDropped_ = new char[numberRows_];
134  memset(rowsDropped_, 0, numberRows_);
135  numberRowsDropped_ = 0;
136  model_ = model;
137  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
138  // Space for starts
139  choleskyStart_ = new CoinBigIndex[numberRows_ + 1];
140  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
141  const int *columnLength = model_->clpMatrix()->getVectorLengths();
142  const int *row = model_->clpMatrix()->getIndices();
143  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
144  const int *rowLength = rowCopy_->getVectorLengths();
145  const int *column = rowCopy_->getIndices();
146  // We need two arrays for counts
147  int *which = new int[numberRows_];
148  int *used = new int[numberRows_ + 1];
149  CoinZeroN(used, numberRows_);
150  int iRow;
151  sizeFactor_ = 0;
152  for (iRow = 0; iRow < numberRows_; iRow++) {
153    int number = 1;
154    // make sure diagonal exists
155    which[0] = iRow;
156    used[iRow] = 1;
157    if (!rowsDropped_[iRow]) {
158      CoinBigIndex startRow = rowStart[iRow];
159      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
160      for (CoinBigIndex k = startRow; k < endRow; k++) {
161        int iColumn = column[k];
162        CoinBigIndex start = columnStart[iColumn];
163        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
164        for (CoinBigIndex j = start; j < end; j++) {
165          int jRow = row[j];
166          if (jRow >= iRow && !rowsDropped_[jRow]) {
167            if (!used[jRow]) {
168              used[jRow] = 1;
169              which[number++] = jRow;
170            }
171          }
172        }
173      }
174      sizeFactor_ += number;
175      int j;
176      for (j = 0; j < number; j++)
177        used[which[j]] = 0;
178    }
179  }
180  delete[] which;
181  // Now we have size - create arrays and fill in
182  try {
183    choleskyRow_ = new int[sizeFactor_];
184  } catch (...) {
185    // no memory
186    delete[] choleskyStart_;
187    choleskyStart_ = NULL;
188    return -1;
189  }
190  try {
191    sparseFactor_ = new double[sizeFactor_];
192  } catch (...) {
193    // no memory
194    delete[] choleskyRow_;
195    choleskyRow_ = NULL;
196    delete[] choleskyStart_;
197    choleskyStart_ = NULL;
198    return -1;
199  }
200
201  sizeFactor_ = 0;
202  which = choleskyRow_;
203  for (iRow = 0; iRow < numberRows_; iRow++) {
204    int number = 1;
205    // make sure diagonal exists
206    which[0] = iRow;
207    used[iRow] = 1;
208    choleskyStart_[iRow] = sizeFactor_;
209    if (!rowsDropped_[iRow]) {
210      CoinBigIndex startRow = rowStart[iRow];
211      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
212      for (CoinBigIndex k = startRow; k < endRow; k++) {
213        int iColumn = column[k];
214        CoinBigIndex start = columnStart[iColumn];
215        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
216        for (CoinBigIndex j = start; j < end; j++) {
217          int jRow = row[j];
218          if (jRow >= iRow && !rowsDropped_[jRow]) {
219            if (!used[jRow]) {
220              used[jRow] = 1;
221              which[number++] = jRow;
222            }
223          }
224        }
225      }
226      sizeFactor_ += number;
227      int j;
228      for (j = 0; j < number; j++)
229        used[which[j]] = 0;
230      // Sort
231      std::sort(which, which + number);
232      // move which on
233      which += number;
234    }
235  }
236  choleskyStart_[numberRows_] = sizeFactor_;
237  delete[] used;
238  permuteInverse_ = new int[numberRows_];
239  permute_ = new int[numberRows_];
240  cholmod_sparse A;
241  A.nrow = numberRows_;
242  A.ncol = numberRows_;
243  A.nzmax = choleskyStart_[numberRows_];
244  A.p = choleskyStart_;
245  A.i = choleskyRow_;
246  A.x = NULL;
247  A.stype = -1;
248  A.itype = CHOLMOD_INT;
249  A.xtype = CHOLMOD_PATTERN;
250  A.dtype = CHOLMOD_DOUBLE;
251  A.sorted = 1;
252  A.packed = 1;
253  c_->nmethods = 9;
254  c_->postorder = true;
255  //c_->dbound=1.0e-20;
256  L_ = cholmod_analyze(&A, c_);
257  if (c_->status) {
258    COIN_DETAIL_PRINT(std::cout << "CHOLMOD ordering failed" << std::endl);
259    return 1;
260  } else {
261    COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", c_->lnz, c_->fl));
262  }
263  for (iRow = 0; iRow < numberRows_; iRow++) {
264    permuteInverse_[iRow] = iRow;
265    permute_[iRow] = iRow;
266  }
267  return 0;
268}
269#endif
270
271/* Does Symbolic factorization given permutation.
272   This is called immediately after order.  If user provides this then
273   user must provide factorize and solve.  Otherwise the default factorization is used
274   returns non-zero if not enough memory */
275int ClpCholeskyUfl::symbolic()
276{
277#ifdef COIN_HAS_CHOLMOD
278  return 0;
279#else
280  return ClpCholeskyBase::symbolic();
281#endif
282}
283
284#ifdef COIN_HAS_CHOLMOD
285/* Factorize - filling in rowsDropped and returning number dropped */
286int ClpCholeskyUfl::factorize(const double *diagonal, int *rowsDropped)
287{
288  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
289  const int *columnLength = model_->clpMatrix()->getVectorLengths();
290  const int *row = model_->clpMatrix()->getIndices();
291  const double *element = model_->clpMatrix()->getElements();
292  const CoinBigIndex *rowStart = rowCopy_->getVectorStarts();
293  const int *rowLength = rowCopy_->getVectorLengths();
294  const int *column = rowCopy_->getIndices();
295  const double *elementByRow = rowCopy_->getElements();
296  int numberColumns = model_->clpMatrix()->getNumCols();
297  int iRow;
298  double *work = new double[numberRows_];
299  CoinZeroN(work, numberRows_);
300  const double *diagonalSlack = diagonal + numberColumns;
301  int newDropped = 0;
302  double largest;
303  //double smallest;
304  //perturbation
305  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
306  perturbation = 0.0;
307  perturbation = perturbation * perturbation;
308  if (perturbation > 1.0) {
309#ifdef COIN_DEVELOP
310    //if (model_->model()->logLevel()&4)
311    std::cout << "large perturbation " << perturbation << std::endl;
312#endif
313    perturbation = sqrt(perturbation);
314    ;
315    perturbation = 1.0;
316  }
317  double delta2 = model_->delta(); // add delta*delta to diagonal
318  delta2 *= delta2;
319  for (iRow = 0; iRow < numberRows_; iRow++) {
320    double *put = sparseFactor_ + choleskyStart_[iRow];
321    int *which = choleskyRow_ + choleskyStart_[iRow];
322    int number = choleskyStart_[iRow + 1] - choleskyStart_[iRow];
323    if (!rowLength[iRow])
324      rowsDropped_[iRow] = 1;
325    if (!rowsDropped_[iRow]) {
326      CoinBigIndex startRow = rowStart[iRow];
327      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
328      work[iRow] = diagonalSlack[iRow] + delta2;
329      for (CoinBigIndex k = startRow; k < endRow; k++) {
330        int iColumn = column[k];
331        if (!whichDense_ || !whichDense_[iColumn]) {
332          CoinBigIndex start = columnStart[iColumn];
333          CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
334          double multiplier = diagonal[iColumn] * elementByRow[k];
335          for (CoinBigIndex j = start; j < end; j++) {
336            int jRow = row[j];
337            if (jRow >= iRow && !rowsDropped_[jRow]) {
338              double value = element[j] * multiplier;
339              work[jRow] += value;
340            }
341          }
342        }
343      }
344      int j;
345      for (j = 0; j < number; j++) {
346        int jRow = which[j];
347        put[j] = work[jRow];
348        work[jRow] = 0.0;
349      }
350    } else {
351      // dropped
352      int j;
353      for (j = 1; j < number; j++) {
354        put[j] = 0.0;
355      }
356      put[0] = 1.0;
357    }
358  }
359  //check sizes
360  double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
361  largest2 *= 1.0e-20;
362  largest = CoinMin(largest2, 1.0e-11);
363  int numberDroppedBefore = 0;
364  for (iRow = 0; iRow < numberRows_; iRow++) {
365    int dropped = rowsDropped_[iRow];
366    // Move to int array
367    rowsDropped[iRow] = dropped;
368    if (!dropped) {
369      CoinBigIndex start = choleskyStart_[iRow];
370      double diagonal = sparseFactor_[start];
371      if (diagonal > largest2) {
372        sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
373      } else {
374        sparseFactor_[start] = CoinMax(diagonal, 1.0e-10);
375        rowsDropped[iRow] = 2;
376        numberDroppedBefore++;
377      }
378    }
379  }
380  delete[] work;
381  cholmod_sparse A;
382  A.nrow = numberRows_;
383  A.ncol = numberRows_;
384  A.nzmax = choleskyStart_[numberRows_];
385  A.p = choleskyStart_;
386  A.i = choleskyRow_;
387  A.x = sparseFactor_;
388  A.stype = -1;
389  A.itype = CHOLMOD_INT;
390  A.xtype = CHOLMOD_REAL;
391  A.dtype = CHOLMOD_DOUBLE;
392  A.sorted = 1;
393  A.packed = 1;
394  cholmod_factorize(&A, L_, c_); /* factorize */
395  choleskyCondition_ = 1.0;
396  bool cleanCholesky;
397  if (model_->numberIterations() < 2000)
398    cleanCholesky = true;
399  else
400    cleanCholesky = false;
401  if (cleanCholesky) {
402    //drop fresh makes some formADAT easier
403    //int oldDropped=numberRowsDropped_;
404    if (newDropped || numberRowsDropped_) {
405      //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
406      //  newDropped<<" dropped)";
407      //if (newDropped>oldDropped)
408      //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
409      //std::cout<<std::endl;
410      newDropped = 0;
411      for (int i = 0; i < numberRows_; i++) {
412        int dropped = rowsDropped[i];
413        rowsDropped_[i] = (char)dropped;
414        if (dropped == 2) {
415          //dropped this time
416          rowsDropped[newDropped++] = i;
417          rowsDropped_[i] = 0;
418        }
419      }
420      numberRowsDropped_ = newDropped;
421      newDropped = -(2 + newDropped);
422    }
423  } else {
424    if (newDropped) {
425      newDropped = 0;
426      for (int i = 0; i < numberRows_; i++) {
427        int dropped = rowsDropped[i];
428        rowsDropped_[i] = (char)dropped;
429        if (dropped == 2) {
430          //dropped this time
431          rowsDropped[newDropped++] = i;
432          rowsDropped_[i] = 1;
433        }
434      }
435    }
436    numberRowsDropped_ += newDropped;
437    if (numberRowsDropped_ && 0) {
438      std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)";
439      if (newDropped) {
440        std::cout << " ( " << newDropped << " dropped this time)";
441      }
442      std::cout << std::endl;
443    }
444  }
445  status_ = 0;
446  return newDropped;
447}
448#else
449/* Factorize - filling in rowsDropped and returning number dropped */
450int ClpCholeskyUfl::factorize(const double *diagonal, int *rowsDropped)
451{
452  return ClpCholeskyBase::factorize(diagonal, rowsDropped);
453}
454#endif
455
456#ifdef COIN_HAS_CHOLMOD
457/* Uses factorization to solve. */
458void ClpCholeskyUfl::solve(double *region)
459{
460  cholmod_dense *x, *b;
461  b = cholmod_allocate_dense(numberRows_, 1, numberRows_, CHOLMOD_REAL, c_);
462  CoinMemcpyN(region, numberRows_, (double *)b->x);
463  x = cholmod_solve(CHOLMOD_A, L_, b, c_);
464  CoinMemcpyN((double *)x->x, numberRows_, region);
465  cholmod_free_dense(&x, c_);
466  cholmod_free_dense(&b, c_);
467}
468#else
469void ClpCholeskyUfl::solve(double *region)
470{
471  ClpCholeskyBase::solve(region);
472}
473#endif
474
475/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
476*/
Note: See TracBrowser for help on using the repository browser.