source: trunk/Clp/src/ClpCholeskyTaucs.cpp @ 2470

Last change on this file since 2470 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.6 KB
Line 
1/* $Id: ClpCholeskyTaucs.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2#ifdef TAUCS_BARRIER
3// Copyright (C) 2004, International Business Machines
4// Corporation and others.  All Rights Reserved.
5// This code is licensed under the terms of the Eclipse Public License (EPL).
6
7#include "CoinPragma.hpp"
8#include "CoinHelperFunctions.hpp"
9#include "ClpHelperFunctions.hpp"
10
11#include "ClpInterior.hpp"
12#include "ClpCholeskyTaucs.hpp"
13#include "ClpMessage.hpp"
14
15//#############################################################################
16// Constructors / Destructor / Assignment
17//#############################################################################
18
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
22ClpCholeskyTaucs::ClpCholeskyTaucs()
23  : ClpCholeskyBase()
24  , matrix_(NULL)
25  , factorization_(NULL)
26  , sparseFactorT_(NULL)
27  , choleskyStartT_(NULL)
28  , choleskyRowT_(NULL)
29  , sizeFactorT_(0)
30  , rowCopyT_(NULL)
31{
32  type_ = 13;
33}
34
35//-------------------------------------------------------------------
36// Copy constructor
37//-------------------------------------------------------------------
38ClpCholeskyTaucs::ClpCholeskyTaucs(const ClpCholeskyTaucs &rhs)
39  : ClpCholeskyBase(rhs)
40{
41  type_ = rhs.type_;
42  // For Taucs stuff is done by malloc
43  matrix_ = rhs.matrix_;
44  sizeFactorT_ = rhs.sizeFactorT_;
45  if (matrix_) {
46    choleskyStartT_ = (int *)malloc((numberRows_ + 1) * sizeof(int));
47    CoinMemcpyN(rhs.choleskyStartT_, (numberRows_ + 1), choleskyStartT_);
48    choleskyRowT_ = (int *)malloc(sizeFactorT_ * sizeof(int));
49    CoinMemcpyN(rhs.choleskyRowT_, sizeFactorT_, choleskyRowT_);
50    sparseFactorT_ = (double *)malloc(sizeFactorT_ * sizeof(double));
51    CoinMemcpyN(rhs.sparseFactorT_, sizeFactorT_, sparseFactorT_);
52    matrix_->colptr = choleskyStartT_;
53    matrix_->rowind = choleskyRowT_;
54    matrix_->values.d = sparseFactorT_;
55  } else {
56    sparseFactorT_ = NULL;
57    choleskyStartT_ = NULL;
58    choleskyRowT_ = NULL;
59  }
60  factorization_ = NULL,
61  rowCopyT_ = rhs.rowCopyT_->clone();
62}
63
64//-------------------------------------------------------------------
65// Destructor
66//-------------------------------------------------------------------
67ClpCholeskyTaucs::~ClpCholeskyTaucs()
68{
69  taucs_ccs_free(matrix_);
70  if (factorization_)
71    taucs_supernodal_factor_free(factorization_);
72  delete rowCopyT_;
73}
74
75//----------------------------------------------------------------
76// Assignment operator
77//-------------------------------------------------------------------
78ClpCholeskyTaucs &
79ClpCholeskyTaucs::operator=(const ClpCholeskyTaucs &rhs)
80{
81  if (this != &rhs) {
82    ClpCholeskyBase::operator=(rhs);
83    taucs_ccs_free(matrix_);
84    if (factorization_)
85      taucs_supernodal_factor_free(factorization_);
86    factorization_ = NULL;
87    sizeFactorT_ = rhs.sizeFactorT_;
88    matrix_ = rhs.matrix_;
89    if (matrix_) {
90      choleskyStartT_ = (int *)malloc((numberRows_ + 1) * sizeof(int));
91      CoinMemcpyN(rhs.choleskyStartT_, (numberRows_ + 1), choleskyStartT_);
92      choleskyRowT_ = (int *)malloc(sizeFactorT_ * sizeof(int));
93      CoinMemcpyN(rhs.choleskyRowT_, sizeFactorT_, choleskyRowT_);
94      sparseFactorT_ = (double *)malloc(sizeFactorT_ * sizeof(double));
95      CoinMemcpyN(rhs.sparseFactorT_, sizeFactorT_, sparseFactorT_);
96      matrix_->colptr = choleskyStartT_;
97      matrix_->rowind = choleskyRowT_;
98      matrix_->values.d = sparseFactorT_;
99    } else {
100      sparseFactorT_ = NULL;
101      choleskyStartT_ = NULL;
102      choleskyRowT_ = NULL;
103    }
104    delete rowCopyT_;
105    rowCopyT_ = rhs.rowCopyT_->clone();
106  }
107  return *this;
108}
109//-------------------------------------------------------------------
110// Clone
111//-------------------------------------------------------------------
112ClpCholeskyBase *ClpCholeskyTaucs::clone() const
113{
114  return new ClpCholeskyTaucs(*this);
115}
116/* Orders rows and saves pointer to matrix.and model */
117int ClpCholeskyTaucs::order(ClpInterior *model)
118{
119  numberRows_ = model->numberRows();
120  rowsDropped_ = new char[numberRows_];
121  memset(rowsDropped_, 0, numberRows_);
122  numberRowsDropped_ = 0;
123  model_ = model;
124  rowCopyT_ = model->clpMatrix()->reverseOrderedCopy();
125  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
126  const int *columnLength = model_->clpMatrix()->getVectorLengths();
127  const int *row = model_->clpMatrix()->getIndices();
128  const CoinBigIndex *rowStart = rowCopyT_->getVectorStarts();
129  const int *rowLength = rowCopyT_->getVectorLengths();
130  const int *column = rowCopyT_->getIndices();
131  // We need two arrays for counts
132  int *which = new int[numberRows_];
133  int *used = new int[numberRows_];
134  CoinZeroN(used, numberRows_);
135  int iRow;
136  sizeFactorT_ = 0;
137  for (iRow = 0; iRow < numberRows_; iRow++) {
138    int number = 1;
139    // make sure diagonal exists
140    which[0] = iRow;
141    used[iRow] = 1;
142    if (!rowsDropped_[iRow]) {
143      CoinBigIndex startRow = rowStart[iRow];
144      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
145      for (CoinBigIndex k = startRow; k < endRow; k++) {
146        int iColumn = column[k];
147        CoinBigIndex start = columnStart[iColumn];
148        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
149        for (CoinBigIndex j = start; j < end; j++) {
150          int jRow = row[j];
151          if (jRow >= iRow && !rowsDropped_[jRow]) {
152            if (!used[jRow]) {
153              used[jRow] = 1;
154              which[number++] = jRow;
155            }
156          }
157        }
158      }
159      sizeFactorT_ += number;
160      int j;
161      for (j = 0; j < number; j++)
162        used[which[j]] = 0;
163    }
164  }
165  delete[] which;
166  // Now we have size - create arrays and fill in
167  matrix_ = taucs_ccs_create(numberRows_, numberRows_, sizeFactorT_,
168    TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
169  if (!matrix_)
170    return 1;
171  // Space for starts
172  choleskyStartT_ = matrix_->colptr;
173  choleskyRowT_ = matrix_->rowind;
174  sparseFactorT_ = matrix_->values.d;
175  sizeFactorT_ = 0;
176  which = choleskyRowT_;
177  for (iRow = 0; iRow < numberRows_; iRow++) {
178    int number = 1;
179    // make sure diagonal exists
180    which[0] = iRow;
181    used[iRow] = 1;
182    choleskyStartT_[iRow] = sizeFactorT_;
183    if (!rowsDropped_[iRow]) {
184      CoinBigIndex startRow = rowStart[iRow];
185      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
186      for (CoinBigIndex k = startRow; k < endRow; k++) {
187        int iColumn = column[k];
188        CoinBigIndex start = columnStart[iColumn];
189        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
190        for (CoinBigIndex j = start; j < end; j++) {
191          int jRow = row[j];
192          if (jRow >= iRow && !rowsDropped_[jRow]) {
193            if (!used[jRow]) {
194              used[jRow] = 1;
195              which[number++] = jRow;
196            }
197          }
198        }
199      }
200      sizeFactorT_ += number;
201      int j;
202      for (j = 0; j < number; j++)
203        used[which[j]] = 0;
204      // Sort
205      std::sort(which, which + number);
206      // move which on
207      which += number;
208    }
209  }
210  choleskyStartT_[numberRows_] = sizeFactorT_;
211  delete[] used;
212  permuteInverse_ = new int[numberRows_];
213  permute_ = new int[numberRows_];
214  int *perm, *invp;
215  // There seem to be bugs in ordering if model too small
216  if (numberRows_ > 10)
217    taucs_ccs_order(matrix_, &perm, &invp, (const char *)"genmmd");
218  else
219    taucs_ccs_order(matrix_, &perm, &invp, (const char *)"identity");
220  CoinMemcpyN(perm, numberRows_, permuteInverse_);
221  free(perm);
222  CoinMemcpyN(invp, numberRows_, permute_);
223  free(invp);
224  // need to permute
225  taucs_ccs_matrix *permuted = taucs_ccs_permute_symmetrically(matrix_, permuteInverse_, permute_);
226  // symbolic
227  factorization_ = taucs_ccs_factor_llt_symbolic(permuted);
228  taucs_ccs_free(permuted);
229  return factorization_ ? 0 : 1;
230}
231/* Does Symbolic factorization given permutation.
232   This is called immediately after order.  If user provides this then
233   user must provide factorize and solve.  Otherwise the default factorization is used
234   returns non-zero if not enough memory */
235int ClpCholeskyTaucs::symbolic()
236{
237  return 0;
238}
239/* Factorize - filling in rowsDropped and returning number dropped */
240int ClpCholeskyTaucs::factorize(const double *diagonal, int *rowsDropped)
241{
242  const CoinBigIndex *columnStart = model_->clpMatrix()->getVectorStarts();
243  const int *columnLength = model_->clpMatrix()->getVectorLengths();
244  const int *row = model_->clpMatrix()->getIndices();
245  const double *element = model_->clpMatrix()->getElements();
246  const CoinBigIndex *rowStart = rowCopyT_->getVectorStarts();
247  const int *rowLength = rowCopyT_->getVectorLengths();
248  const int *column = rowCopyT_->getIndices();
249  const double *elementByRow = rowCopyT_->getElements();
250  int numberColumns = model_->clpMatrix()->getNumCols();
251  int iRow;
252  double *work = new double[numberRows_];
253  CoinZeroN(work, numberRows_);
254  const double *diagonalSlack = diagonal + numberColumns;
255  int newDropped = 0;
256  double largest;
257  //perturbation
258  double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
259  perturbation = perturbation * perturbation;
260  if (perturbation > 1.0) {
261    //if (model_->model()->logLevel()&4)
262    std::cout << "large perturbation " << perturbation << std::endl;
263    perturbation = sqrt(perturbation);
264    ;
265    perturbation = 1.0;
266  }
267  for (iRow = 0; iRow < numberRows_; iRow++) {
268    double *put = sparseFactorT_ + choleskyStartT_[iRow];
269    int *which = choleskyRowT_ + choleskyStartT_[iRow];
270    int number = choleskyStartT_[iRow + 1] - choleskyStartT_[iRow];
271    if (!rowLength[iRow])
272      rowsDropped_[iRow] = 1;
273    if (!rowsDropped_[iRow]) {
274      CoinBigIndex startRow = rowStart[iRow];
275      CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
276      work[iRow] = diagonalSlack[iRow];
277      for (CoinBigIndex k = startRow; k < endRow; k++) {
278        int iColumn = column[k];
279        CoinBigIndex start = columnStart[iColumn];
280        CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
281        double multiplier = diagonal[iColumn] * elementByRow[k];
282        for (CoinBigIndex j = start; j < end; j++) {
283          int jRow = row[j];
284          if (jRow >= iRow && !rowsDropped_[jRow]) {
285            double value = element[j] * multiplier;
286            work[jRow] += value;
287          }
288        }
289      }
290      int j;
291      for (j = 0; j < number; j++) {
292        int jRow = which[j];
293        put[j] = work[jRow];
294        work[jRow] = 0.0;
295      }
296    } else {
297      // dropped
298      int j;
299      for (j = 1; j < number; j++) {
300        put[j] = 0.0;
301      }
302      put[0] = 1.0;
303    }
304  }
305  //check sizes
306  double largest2 = maximumAbsElement(sparseFactorT_, sizeFactorT_);
307  largest2 *= 1.0e-19;
308  largest = CoinMin(largest2, 1.0e-11);
309  int numberDroppedBefore = 0;
310  for (iRow = 0; iRow < numberRows_; iRow++) {
311    int dropped = rowsDropped_[iRow];
312    // Move to int array
313    rowsDropped[iRow] = dropped;
314    if (!dropped) {
315      CoinBigIndex start = choleskyStartT_[iRow];
316      double diagonal = sparseFactorT_[start];
317      if (diagonal > largest2) {
318        sparseFactorT_[start] = diagonal + perturbation;
319      } else {
320        sparseFactorT_[start] = diagonal + perturbation;
321        rowsDropped[iRow] = 2;
322        numberDroppedBefore++;
323      }
324    }
325  }
326  taucs_supernodal_factor_free_numeric(factorization_);
327  // need to permute
328  taucs_ccs_matrix *permuted = taucs_ccs_permute_symmetrically(matrix_, permuteInverse_, permute_);
329  int rCode = taucs_ccs_factor_llt_numeric(permuted, factorization_);
330  taucs_ccs_free(permuted);
331  if (rCode)
332    printf("return code of %d from factor\n", rCode);
333  delete[] work;
334  choleskyCondition_ = 1.0;
335  bool cleanCholesky;
336  if (model_->numberIterations() < 200)
337    cleanCholesky = true;
338  else
339    cleanCholesky = false;
340  /*
341       How do I find out where 1.0e100's are in cholesky?
342     */
343  if (cleanCholesky) {
344    //drop fresh makes some formADAT easier
345    int oldDropped = numberRowsDropped_;
346    if (newDropped || numberRowsDropped_) {
347      std::cout << "Rank " << numberRows_ - newDropped << " ( " << newDropped << " dropped)";
348      if (newDropped > oldDropped)
349        std::cout << " ( " << newDropped - oldDropped << " dropped this time)";
350      std::cout << std::endl;
351      newDropped = 0;
352      for (int i = 0; i < numberRows_; i++) {
353        char dropped = rowsDropped[i];
354        rowsDropped_[i] = dropped;
355        if (dropped == 2) {
356          //dropped this time
357          rowsDropped[newDropped++] = i;
358          rowsDropped_[i] = 0;
359        }
360      }
361      numberRowsDropped_ = newDropped;
362      newDropped = -(1 + newDropped);
363    }
364  } else {
365    if (newDropped) {
366      newDropped = 0;
367      for (int i = 0; i < numberRows_; i++) {
368        char dropped = rowsDropped[i];
369        int oldDropped = rowsDropped_[i];
370        ;
371        rowsDropped_[i] = dropped;
372        if (dropped == 2) {
373          assert(!oldDropped);
374          //dropped this time
375          rowsDropped[newDropped++] = i;
376          rowsDropped_[i] = 1;
377        }
378      }
379    }
380    numberRowsDropped_ += newDropped;
381    if (numberRowsDropped_) {
382      std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)";
383      if (newDropped) {
384        std::cout << " ( " << newDropped << " dropped this time)";
385      }
386      std::cout << std::endl;
387    }
388  }
389  status_ = 0;
390  return newDropped;
391}
392/* Uses factorization to solve. */
393void ClpCholeskyTaucs::solve(double *region)
394{
395  double *in = new double[numberRows_];
396  double *out = new double[numberRows_];
397  taucs_vec_permute(numberRows_, TAUCS_DOUBLE, region, in, permuteInverse_);
398  int rCode = taucs_supernodal_solve_llt(factorization_, out, in);
399  if (rCode)
400    printf("return code of %d from solve\n", rCode);
401  taucs_vec_permute(numberRows_, TAUCS_DOUBLE, out, region, permute_);
402  delete[] out;
403  delete[] in;
404}
405#endif
406
407/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
408*/
Note: See TracBrowser for help on using the repository browser.