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
412