1/* $Id: CoinPresolveMatrix.cpp 1448 2011-06-19 15:34:41Z 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#include <stdio.h>
7
8#include <cassert>
9#include <iostream>
10
11#include "CoinHelperFunctions.hpp"
12#include "CoinPresolveMatrix.hpp"
13#include "CoinTime.hpp"
14
15/*! \file
16
17 This file contains methods for CoinPresolveMatrix, the object used during
18 presolve transformations.
19*/
20
21/*
22 Constructor and destructor for CoinPresolveMatrix.
23*/
24
25/*
26 CoinPresolveMatrix constructor
27
28 The constructor does very little, for much the same reasons that the
29 CoinPrePostsolveMatrix constructor does little. Might as well wait until we
30 load a matrix.
31
32 In general, for presolve the allocated size can be equal to the size of the
33 constraint matrix before presolve transforms are applied. (Presolve
34 transforms are assumed to reduce the size of the constraint system.) But we
35 need to keep the *_alloc parameters for compatibility with
36 CoinPrePostsolveMatrix.
37*/
38
39CoinPresolveMatrix::CoinPresolveMatrix
40 (int ncols_alloc, int nrows_alloc, CoinBigIndex nelems_alloc)
41
42 : CoinPrePostsolveMatrix(ncols_alloc,nrows_alloc,nelems_alloc),
43
44 clink_(0),
45 rlink_(0),
46
47 dobias_(0.0),
48 mrstrt_(0),
49 hinrow_(0),
50 rowels_(0),
51 hcol_(0),
52
53 integerType_(0),
54 anyInteger_(false),
55 tuning_(false),
56 startTime_(0.0),
57 feasibilityTolerance_(0.0),
58 status_(-1),
59 pass_(0),
60 maxSubstLevel_(3),
61 colChanged_(0),
62 colsToDo_(0),
63 numberColsToDo_(0),
64 nextColsToDo_(0),
65 numberNextColsToDo_(0),
66
67 rowChanged_(0),
68 rowsToDo_(0),
69 numberRowsToDo_(0),
70 nextRowsToDo_(0),
71 numberNextRowsToDo_(0),
72 presolveOptions_(0),
73 anyProhibited_(false),
74 usefulRowInt_(NULL),
75 usefulRowDouble_(NULL),
76 usefulColumnInt_(NULL),
77 usefulColumnDouble_(NULL),
78 randomNumber_(NULL),
79 infiniteUp_(NULL),
80 sumUp_(NULL),
81 infiniteDown_(NULL),
82 sumDown_(NULL)
83
84{ /* nothing to do here */
85
86 return ; }
87
88/*
89 CoinPresolveMatrix destructor.
90*/
91
92CoinPresolveMatrix::~CoinPresolveMatrix()
93
94{ delete[] clink_ ;
95 delete[] rlink_ ;
96
97 delete[] mrstrt_ ;
98 delete[] hinrow_ ;
99 delete[] rowels_ ;
100 delete[] hcol_ ;
101
102 delete[] integerType_ ;
103 delete[] rowChanged_ ;
104 delete[] rowsToDo_ ;
105 delete[] nextRowsToDo_ ;
106 delete[] colChanged_ ;
107 delete[] colsToDo_ ;
108 delete[] nextColsToDo_ ;
109 delete[] usefulRowInt_;
110 delete[] usefulRowDouble_;
111 delete[] usefulColumnInt_;
112 delete[] usefulColumnDouble_;
113 delete[] randomNumber_;
114 delete[] infiniteUp_;
115 delete[] sumUp_;
116 delete[] infiniteDown_;
117 delete[] sumDown_;
118
119 return ; }
120
121
122
123/*
124 This routine loads a CoinPackedMatrix and proceeds to do the bulk of the
125 initialisation for the PrePostsolve and Presolve objects.
126
127 In the CoinPrePostsolveMatrix portion of the object, it initialises the
128 column-major packed matrix representation and the arrays that track the
129 motion of original columns and rows.
130
131 In the CoinPresolveMatrix portion of the object, it initialises the
132 row-major packed matrix representation, the arrays that assist in matrix
133 storage management, and the arrays that track the rows and columns to be
134 processed.
135
136 Arrays are allocated to the requested size (ncols0_, nrow0_, nelems0_).
137
138 The source matrix must be column ordered; it does not need to be gap-free.
139 Bulk storage in the column-major (hrow_, colels_) and row-major (hcol_,
140 rowels_) matrices is allocated at twice the required size so that we can
141 expand columns and rows as needed. This is almost certainly grossly
142 oversize, but (1) it's efficient, and (2) the utility routines which
143 compact the bulk storage areas have no provision to reallocate.
144*/
145
146void CoinPresolveMatrix::setMatrix (const CoinPackedMatrix *mtx)
147
148{
149/*
150 Check to make sure the matrix will fit and is column ordered.
151*/
152 if (mtx->isColOrdered() == false)
153 { throw CoinError("source matrix must be column ordered",
154 "setMatrix","CoinPrePostsolveMatrix") ; }
155
156 int numCols = mtx->getNumCols() ;
157 if (numCols > ncols0_)
158 { throw CoinError("source matrix exceeds allocated capacity",
159 "setMatrix","CoinPrePostsolveMatrix") ; }
160/*
161 Acquire the actual size, but allocate the matrix storage to the
162 requested capacity. The column-major rep is part of the PrePostsolve
163 object, the row-major rep belongs to the Presolve object.
164*/
165 ncols_ = numCols ;
166 nrows_ = mtx->getNumRows() ;
167 nelems_ = mtx->getNumElements() ;
168 bulk0_ = static_cast<CoinBigIndex> (bulkRatio_*nelems0_) ;
169
170 if (mcstrt_ == 0) mcstrt_ = new CoinBigIndex [ncols0_+1] ;
171 if (hincol_ == 0) hincol_ = new int [ncols0_+1] ;
172 if (hrow_ == 0) hrow_ = new int [bulk0_] ;
173 if (colels_ == 0) colels_ = new double [bulk0_] ;
174
175 if (mrstrt_ == 0) mrstrt_ = new CoinBigIndex [nrows0_+1] ;
176 if (hinrow_ == 0) hinrow_ = new int [nrows0_+1] ;
177 if (hcol_ == 0) hcol_ = new int [bulk0_] ;
178 if (rowels_ == 0) rowels_ = new double [bulk0_] ;
179/*
180 Grab the corresponding vectors from the source matrix.
181*/
182 const CoinBigIndex *src_mcstrt = mtx->getVectorStarts() ;
183 const int *src_hincol = mtx->getVectorLengths() ;
184 const double *src_colels = mtx->getElements() ;
185 const int *src_hrow = mtx->getIndices() ;
186/*
187 Bulk copy the column starts and lengths.
188*/
189 CoinMemcpyN(src_mcstrt,mtx->getSizeVectorStarts(),mcstrt_) ;
190 CoinMemcpyN(src_hincol,mtx->getSizeVectorLengths(),hincol_) ;
191/*
192 Copy the coefficients column by column in case there are gaps between
193 the columns in the bulk storage area. The assert is just in case the
194 gaps are *really* big.
195*/
196 assert(src_mcstrt[ncols_] <= bulk0_) ;
197 int j;
198 for ( j = 0 ; j < numCols ; j++)
199 { int lenj = src_hincol[j] ;
200 CoinBigIndex offset = mcstrt_[j] ;
201 CoinMemcpyN(src_colels+offset,lenj,colels_+offset) ;
202 CoinMemcpyN(src_hrow+offset,lenj,hrow_+offset) ; }
203/*
204 Now make a row-major copy. Start by counting the number of coefficients in
205 each row; we can do this directly in hinrow. Given the number of
206 coefficients in a row, we know how to lay out the bulk storage area.
207*/
208 CoinZeroN(hinrow_,nrows0_+1) ;
209 for ( j = 0 ; j < ncols_ ; j++)
210 { int *rowIndices = hrow_+mcstrt_[j] ;
211 int lenj = hincol_[j] ;
212 for (int k = 0 ; k < lenj ; k++)
213 { int i = rowIndices[k] ;
214 hinrow_[i]++ ; } }
215/*
216 Initialize mrstrt[i] to the start of row i+1. As we drop each coefficient
217 and column index into the bulk storage arrays, we'll decrement and store.
218 When we're done, mrstrt[i] will point to the start of row i, as it should.
219*/
220 int totalCoeffs = 0 ;
221 int i;
222 for ( i = 0 ; i < nrows_ ; i++)
223 { totalCoeffs += hinrow_[i] ;
224 mrstrt_[i] = totalCoeffs ; }
225 mrstrt_[nrows_] = totalCoeffs ;
226 for ( j = ncols_-1 ; j >= 0 ; j--)
227 { int lenj = hincol_[j] ;
228 double *colCoeffs = colels_+mcstrt_[j] ;
229 int *rowIndices = hrow_+mcstrt_[j] ;
230 for (int k = 0 ; k < lenj ; k++)
231 { int ri;
232 ri = rowIndices[k] ;
233 double aij = colCoeffs[k] ;
234 CoinBigIndex l = --mrstrt_[ri] ;
235 rowels_[l] = aij ;
236 hcol_[l] = j ; } }
237/*
238 Now the support structures. The entry for original column j should start
239 out as j; similarly for row i. originalColumn_ and originalRow_ belong to
240 the PrePostsolve object.
241*/
242 if (originalColumn_ == 0) originalColumn_ = new int [ncols0_] ;
243 if (originalRow_ == 0) originalRow_ = new int [nrows0_] ;
244
245 for ( j = 0 ; j < ncols0_ ; j++)
246 originalColumn_[j] = j ;
247 for ( i = 0 ; i < nrows0_ ; i++)
248 originalRow_[i] = i ;
249/*
250 We have help to set up the clink_ and rlink_ vectors (aids for matrix bulk
251 storage management). clink_ and rlink_ belong to the Presolve object. Once
252 this is done, it's safe to set mrstrt_[nrows_] and mcstrt_[ncols_] to the
253 full size of the bulk storage area.
254*/
255 if (clink_ == 0) clink_ = new presolvehlink [ncols0_+1] ;
256 if (rlink_ == 0) rlink_ = new presolvehlink [nrows0_+1] ;
257 presolve_make_memlists(/*mcstrt_,*/hincol_,clink_,ncols_) ;
258 presolve_make_memlists(/*mrstrt_,*/hinrow_,rlink_,nrows_) ;
259 mcstrt_[ncols_] = bulk0_ ;
260 mrstrt_[nrows_] = bulk0_ ;
261/*
262 No rows or columns have been changed just yet. colChanged_ and rowChanged_
263 belong to the Presolve object.
264*/
265 if (colChanged_ == 0) colChanged_ = new unsigned char [ncols0_] ;
266 CoinZeroN(colChanged_,ncols0_) ;
267 if (rowChanged_ == 0) rowChanged_ = new unsigned char [nrows0_] ;
268 CoinZeroN(rowChanged_,nrows0_) ;
269/*
270 Finally, allocate the various *ToDo arrays. These are used to track the rows
271 and columns which should be processed in a given round of presolve
272 transforms. These belong to the Presolve object. Setting number*ToDo to 0
273 is all the initialization that's required here.
274*/
275 rowsToDo_ = new int [nrows0_] ;
276 numberRowsToDo_ = 0 ;
277 nextRowsToDo_ = new int [nrows0_] ;
278 numberNextRowsToDo_ = 0 ;
279 colsToDo_ = new int [ncols0_] ;
280 numberColsToDo_ = 0 ;
281 nextColsToDo_ = new int [ncols0_] ;
282 numberNextColsToDo_ = 0 ;
283 initializeStuff();
284 return ; }
285/* Recompute ups and downs for a row (nonzero if infeasible).
286 If iRow -1 then recompute all */
287int
288CoinPresolveMatrix::recomputeSums(int iRow)
289{
290 double * columnLower = clo_;
291 double * columnUpper = cup_;
292
293 const double *element = rowels_;
294 const int *column = hcol_;
295 const CoinBigIndex *rowStart = mrstrt_;
296 const int *rowLength = hinrow_;
297 int numberRows = nrows_;
298 int numberColumns = ncols_;
299 //const int *hrow = hrow_;
300 //const CoinBigIndex *mcstrt = mcstrt_;
301 //const int *hincol = hincol_;
302 double *rowLower = rlo_;
303 double *rowUpper = rup_;
304 double large = 1.0e20; // treat bounds > this as infinite
305 int iFirst = (iRow>=0) ? iRow : 0;
306 int iLast = (iRow>=0) ? iRow : numberRows;
307 int infeasible=0;
308 double tolerance = feasibilityTolerance_;
309 for (iRow=iFirst;iRow<iLast;iRow++) {
310 infiniteUp_[iRow]=0;
311 sumUp_[iRow]=0.0;
312 infiniteDown_[iRow]=0;
313 sumDown_[iRow]=0.0;
314 if ((rowLower[iRow]>-large||rowUpper[iRow]<large)&&rowLength[iRow]>0) {
315 int infiniteUpper = 0;
316 int infiniteLower = 0;
317 double maximumUp = 0.0;
318 double maximumDown = 0.0;
319 CoinBigIndex rStart = rowStart[iRow];
320 CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
321 CoinBigIndex j;
322 // Compute possible lower and upper ranges
323
324 for (j = rStart; j < rEnd; ++j) {
325 double value=element[j];
326 int iColumn = column[j];
327 if (value > 0.0) {
328 if (columnUpper[iColumn] < large)
329 maximumUp += columnUpper[iColumn] * value;
330 else
331 ++infiniteUpper;
332 if (columnLower[iColumn] > -large)
333 maximumDown += columnLower[iColumn] * value;
334 else
335 ++infiniteLower;
336 } else if (value<0.0) {
337 if (columnUpper[iColumn] < large)
338 maximumDown += columnUpper[iColumn] * value;
339 else
340 ++infiniteLower;
341 if (columnLower[iColumn] > -large)
342 maximumUp += columnLower[iColumn] * value;
343 else
344 ++infiniteUpper;
345 }
346 }
347#if 0
348 // Build in a margin of error (NO)
349 maximumUp += 1.0e-8*fabs(maximumUp);
350 maximumDown -= 1.0e-8*fabs(maximumDown);
351#endif
352 infiniteUp_[iRow]=infiniteUpper;
353 sumUp_[iRow]=maximumUp;
354 infiniteDown_[iRow]=infiniteLower;
355 sumDown_[iRow]=maximumDown;
356 double maxUp = maximumUp+infiniteUpper*1.0e31;
357 double maxDown = maximumDown-infiniteLower*1.0e31;
358 if (maxUp <= rowUpper[iRow] + tolerance &&
359 maxDown >= rowLower[iRow] - tolerance) {
360 // redundant
361 infiniteUp_[iRow]=numberColumns+1;
362 infiniteDown_[iRow]=numberColumns+1;
363 } else if (maxUp <rowLower[iRow]-tolerance) {
364 // infeasible
365 infeasible++;
366 } else if (maxDown >rowUpper[iRow]+tolerance) {
367 // infeasible
368 infeasible++;
369 }
370 } else {
371 // odd probably redundant
372 infiniteUp_[iRow]=numberColumns+1;
373 infiniteDown_[iRow]=numberColumns+1;
374 if (rowLower[iRow]>0.0||rowUpper[iRow]<0.0) {
375 double tolerance2=10.0*tolerance;
376 if (rowLower[iRow]>0.0&&rowLower[iRow]<tolerance2)
377 rowLower[iRow]=0.0;
378 else
379 infeasible++;
380 if (rowUpper[iRow]<0.0&&rowUpper[iRow]>-tolerance2)
381 rowUpper[iRow]=0.0;
382 else
383 infeasible++;
384 }
385 }
386 }
387 return infeasible;
388}
389// Initialize random numbers etc (nonzero if infeasible)
390int
391CoinPresolveMatrix::initializeStuff()
392{
393 // Allocate useful arrays
394 usefulRowInt_ = new int [3*nrows_];
395 usefulRowDouble_ = new double [nrows_];
396 usefulColumnInt_ = new int [2*ncols_];
397 usefulColumnDouble_ = new double[ncols_];
398 int k=CoinMax(ncols_+1,nrows_+1);
399 randomNumber_ = new double [k];
400 coin_init_random_vec(randomNumber_,k);
401 infiniteUp_ = new int [nrows_];
402 sumUp_ = new double [nrows_];
403 infiniteDown_ = new int [nrows_];
404 sumDown_ = new double [nrows_];
405 // return recomputeSums(-1);
406 return 0;
407}
408// Delete useful arrays
409void
410CoinPresolveMatrix::deleteStuff()
411{
412 delete[] usefulRowInt_;
413 delete[] usefulRowDouble_;
414 delete[] usefulColumnInt_;
415 delete[] usefulColumnDouble_;
416 delete[] randomNumber_;
417 delete[] infiniteUp_;
418 delete[] sumUp_;
419 delete[] infiniteDown_;
420 delete[] sumDown_;
421 usefulRowInt_ = NULL;
422 usefulRowDouble_ = NULL;
423 usefulColumnInt_ = NULL;
424 usefulColumnDouble_ = NULL;
425 randomNumber_ = NULL;
426 infiniteUp_ = NULL;
427 sumUp_ = NULL;
428 infiniteDown_ = NULL;
429 sumDown_ = NULL;
430}
431
432
433/*
434 These functions set integer type information. The first expects an array with
435 an entry for each variable. The second sets all variables to integer or
436 continuous type.
437*/
438
439void CoinPresolveMatrix::setVariableType (const unsigned char *variableType,
440 int lenParam)
441
442{ int len ;
443
444 if (lenParam < 0)
445 { len = ncols_ ; }
446 else
447 if (lenParam > ncols0_)
448 { throw CoinError("length exceeds allocated size",
449 "setIntegerType","CoinPresolveMatrix") ; }
450 else
451 { len = lenParam ; }
452
453 if (integerType_ == 0) integerType_ = new unsigned char [ncols0_] ;
454 CoinCopyN(variableType,len,integerType_) ;
455
456 return ; }
457
458void CoinPresolveMatrix::setVariableType (bool allIntegers, int lenParam)
459
460{ int len ;
461
462 if (lenParam < 0)
463 { len = ncols_ ; }
464 else
465 if (lenParam > ncols0_)
466 { throw CoinError("length exceeds allocated size",
467 "setIntegerType","CoinPresolveMatrix") ; }
468 else
469 { len = lenParam ; }
470
471 if (integerType_ == 0) integerType_ = new unsigned char [ncols0_] ;
472
473 const unsigned char value = 1 ;
474
475 if (allIntegers == true)
476 { CoinFillN(integerType_,len,value) ; }
477 else
478 { CoinZeroN(integerType_,len) ; }
479
480 return ; }
481
482/*
483 The next pair of routines initialises the [row,col]ToDo lists in preparation
484 for a major pass. All except rows/columns marked as prohibited are added to
485 the lists.
486*/
487
488void CoinPresolveMatrix::initColsToDo ()
489/*
490 Initialize the ToDo lists in preparation for a major iteration of
491 preprocessing. First, cut back the ToDo and NextToDo lists to zero entries.
492 Then place all columns not marked prohibited on the ToDo list.
493*/
494
495{ int j ;
496
497 numberNextColsToDo_ = 0 ;
498
499 if (anyProhibited_ == false)
500 { for (j = 0 ; j < ncols_ ; j++)
501 { colsToDo_[j] = j ; }
502 numberColsToDo_ = ncols_ ; }
503 else
504 { numberColsToDo_ = 0 ;
505 for (j = 0 ; j < ncols_ ; j++)
506 if (colProhibited(j) == false)
507 { colsToDo_[numberColsToDo_++] = j ; } }
508
509 return ; }
510
511void CoinPresolveMatrix::initRowsToDo ()
512/*
513 Initialize the ToDo lists in preparation for a major iteration of
514 preprocessing. First, cut back the ToDo and NextToDo lists to zero entries.
515 Then place all rows not marked prohibited on the ToDo list.
516*/
517
518{ int i ;
519
520 numberNextRowsToDo_ = 0 ;
521
522 if (anyProhibited_ == false)
523 { for (i = 0 ; i < nrows_ ; i++)
524 { rowsToDo_[i] = i ; }
525 numberRowsToDo_ = nrows_ ; }
526 else
527 { numberRowsToDo_ = 0 ;
528 for (i = 0 ; i < nrows_ ; i++)
529 if (rowProhibited(i) == false)
530 { rowsToDo_[numberRowsToDo_++] = i ; } }
531
532 return ; }
533
534int CoinPresolveMatrix::stepColsToDo ()
535/*
536 This routine transfers the contents of NextToDo to ToDo, simultaneously
537 resetting the Changed indicator. It returns the number of columns
538 transfered.
539*/
540{ int k ;
541
542 for (k = 0 ; k < numberNextColsToDo_ ; k++)
543 { int j = nextColsToDo_[k] ;
544 unsetColChanged(j) ;
545 colsToDo_[k] = j ; }
546 numberColsToDo_ = numberNextColsToDo_ ;
547 numberNextColsToDo_ = 0 ;
548
549 return (numberColsToDo_) ; }
550
551int CoinPresolveMatrix::stepRowsToDo ()
552/*
553 This routine transfers the contents of NextToDo to ToDo, simultaneously
554 resetting the Changed indicator. It returns the number of columns
555 transfered.
556*/
557{ int k ;
558
559 for (k = 0 ; k < numberNextRowsToDo_ ; k++)
560 { int i = nextRowsToDo_[k] ;
561 unsetRowChanged(i) ;
562 rowsToDo_[k] = i ; }
563 numberRowsToDo_ = numberNextRowsToDo_ ;
564 numberNextRowsToDo_ = 0 ;
565
566 return (numberRowsToDo_) ; }
567// Say we want statistics - also set time
568void
569CoinPresolveMatrix::statistics()
570{
571 tuning_=true;
572 startTime_ = CoinCpuTime();
573}
574#ifdef PRESOLVE_DEBUG
575#include "CoinPresolvePsdebug.cpp"
576#endif
577