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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

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