1/* $Id: CoinPrePostsolveMatrix.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 <cstdio>
7#include <cassert>
8#include <iostream>
9
10#include "CoinHelperFunctions.hpp"
11#include "CoinPresolveMatrix.hpp"
12#ifndef SLIM_CLP
13#include "CoinWarmStartBasis.hpp"
14#endif
15
16/*! \file
17 This file contains methods for CoinPrePostsolveMatrix, the foundation class
18 for CoinPresolveMatrix and CoinPostsolveMatrix.
19*/
20
21/*
22 Constructor and destructor for CoinPrePostsolveMatrix.
23*/
24
25/*
26 CoinPrePostsolveMatrix constructor
27
28 This constructor does next to nothing, because there's no sensible middle
29 ground between next to nothing and a constructor with twenty parameters
30 that all need to be extracted from the constraint system held by an OSI.
31 The alternative, creating a constructor which takes some flavour of OSI as
32 a parameter, seems to me (lh) to be wrong. That knowledge does not belong
33 in the generic COIN support library.
34
35 The philosophy here is to create an empty CoinPrePostsolveMatrix object and
36 then load in the constraint matrix, vectors, and miscellaneous parameters.
37 Some of this will be done from CoinPresolveMatrix or CoinPostsolveMatrix
38 constructors, but in the end most of it should be pushed back to an
39 OSI-specific method. Then the knowledge of how to access the required data
40 in an OSI is pushed back to the individual OSI classes where it belongs.
41
42 Thus, all vector allocation is postponed until load time.
43*/
44
45CoinPrePostsolveMatrix::CoinPrePostsolveMatrix
46 (int ncols_alloc, int nrows_alloc, CoinBigIndex nelems_alloc)
47
48 : ncols_(0),
49 nrows_(0),
50 nelems_(0),
51 ncols0_(ncols_alloc),
52 nrows0_(nrows_alloc),
53 nelems0_(nelems_alloc),
54 bulkRatio_(2.0),
55
56 mcstrt_(0),
57 hincol_(0),
58 hrow_(0),
59 colels_(0),
60
61 cost_(0),
62 originalOffset_(0),
63 clo_(0),
64 cup_(0),
65 rlo_(0),
66 rup_(0),
67
68 originalColumn_(0),
69 originalRow_(0),
70
71 ztolzb_(0.0),
72 ztoldj_(0.0),
73
74 maxmin_(0),
75
76 sol_(0),
77 rowduals_(0),
78 acts_(0),
79 rcosts_(0),
80 colstat_(0),
81 rowstat_(0),
82
83 handler_(0),
84 defaultHandler_(false),
85 messages_()
86
87{ handler_ = new CoinMessageHandler() ;
88 defaultHandler_ = true ;
89 bulk0_ = static_cast<CoinBigIndex> (bulkRatio_*nelems_alloc);
90
91 return ; }
92
93/*
94 CoinPrePostsolveMatrix destructor
95*/
96
97CoinPrePostsolveMatrix::~CoinPrePostsolveMatrix()
98{
99 delete[] sol_ ;
100 delete[] rowduals_ ;
101 delete[] acts_ ;
102 delete[] rcosts_ ;
103
104/*
105 Note that we do NOT delete rowstat_. This is to maintain compatibility with
106 ClpPresolve and OsiPresolve, which allocate a single vector and split it
107 between column and row status.
108*/
109 delete[] colstat_ ;
110
111 delete[] cost_ ;
112 delete[] clo_ ;
113 delete[] cup_ ;
114 delete[] rlo_ ;
115 delete[] rup_ ;
116
117 delete[] mcstrt_ ;
118 delete[] hrow_ ;
119 delete[] colels_ ;
120 delete[] hincol_ ;
121
122 delete[] originalColumn_ ;
123 delete[] originalRow_ ;
124
125 if (defaultHandler_ == true)
126 delete handler_ ;
127}
128
129
130#ifndef SLIM_CLP
131/*
132 Methods to set the miscellaneous parameters: max/min, objective offset, and
133 tolerances.
134*/
135
136void CoinPrePostsolveMatrix::setObjOffset (double offset)
137
138{ originalOffset_ = offset ; }
139
140void CoinPrePostsolveMatrix::setObjSense (double objSense)
141
142{ maxmin_ = objSense ; }
143
144void CoinPrePostsolveMatrix::setPrimalTolerance (double primTol)
145
146{ ztolzb_ = primTol ; }
147
148void CoinPrePostsolveMatrix::setDualTolerance (double dualTol)
149
150{ ztoldj_ = dualTol ; }
151
152
153
154/*
155 Methods to set the various vectors. For all methods, lenParam can be
156 omitted and will default to -1. In that case, the default action is to copy
157 ncols_ or nrows_ entries, as appropriate.
158
159 It is *not* considered an error to specify lenParam = 0! This allows for
160 allocation of vectors in an intially empty system. Note that ncols_ and
161 nrows_ will be 0 before a constraint system is loaded. Be careful what you
162 ask for.
163
164 The vector allocated in the CoinPrePostsolveMatrix will be of size ncols0_
165 or nrows0_, as appropriate.
166*/
167
168void CoinPrePostsolveMatrix::setColLower (const double *colLower, int lenParam)
169
170{ int len ;
171
172 if (lenParam < 0)
173 { len = ncols_ ; }
174 else
175 if (lenParam > ncols0_)
176 { throw CoinError("length exceeds allocated size",
177 "setColLower","CoinPrePostsolveMatrix") ; }
178 else
179 { len = lenParam ; }
180
181 if (clo_ == 0) clo_ = new double[ncols0_] ;
182 CoinMemcpyN(colLower,len,clo_) ;
183
184 return ; }
185
186void CoinPrePostsolveMatrix::setColUpper (const double *colUpper, int lenParam)
187
188{ int len ;
189
190 if (lenParam < 0)
191 { len = ncols_ ; }
192 else
193 if (lenParam > ncols0_)
194 { throw CoinError("length exceeds allocated size",
195 "setColUpper","CoinPrePostsolveMatrix") ; }
196 else
197 { len = lenParam ; }
198
199 if (cup_ == 0) cup_ = new double[ncols0_] ;
200 CoinMemcpyN(colUpper,len,cup_) ;
201
202 return ; }
203
204void CoinPrePostsolveMatrix::setColSolution (const double *colSol,
205 int lenParam)
206
207{ int len ;
208
209 if (lenParam < 0)
210 { len = ncols_ ; }
211 else
212 if (lenParam > ncols0_)
213 { throw CoinError("length exceeds allocated size",
214 "setColSolution","CoinPrePostsolveMatrix") ; }
215 else
216 { len = lenParam ; }
217
218 if (sol_ == 0) sol_ = new double[ncols0_] ;
219 CoinMemcpyN(colSol,len,sol_) ;
220
221 return ; }
222
223void CoinPrePostsolveMatrix::setCost (const double *cost, int lenParam)
224
225{ int len ;
226
227 if (lenParam < 0)
228 { len = ncols_ ; }
229 else
230 if (lenParam > ncols0_)
231 { throw CoinError("length exceeds allocated size",
232 "setCost","CoinPrePostsolveMatrix") ; }
233 else
234 { len = lenParam ; }
235
236 if (cost_ == 0) cost_ = new double[ncols0_] ;
237 CoinMemcpyN(cost,len,cost_) ;
238
239 return ; }
240
241void CoinPrePostsolveMatrix::setReducedCost (const double *redCost,
242 int lenParam)
243
244{ int len ;
245
246 if (lenParam < 0)
247 { len = ncols_ ; }
248 else
249 if (lenParam > ncols0_)
250 { throw CoinError("length exceeds allocated size",
251 "setReducedCost","CoinPrePostsolveMatrix") ; }
252 else
253 { len = lenParam ; }
254
255 if (rcosts_ == 0) rcosts_ = new double[ncols0_] ;
256 CoinMemcpyN(redCost,len,rcosts_) ;
257
258 return ; }
259
260
261void CoinPrePostsolveMatrix::setRowLower (const double *rowLower, int lenParam)
262
263{ int len ;
264
265 if (lenParam < 0)
266 { len = nrows_ ; }
267 else
268 if (lenParam > nrows0_)
269 { throw CoinError("length exceeds allocated size",
270 "setRowLower","CoinPrePostsolveMatrix") ; }
271 else
272 { len = lenParam ; }
273
274 if (rlo_ == 0) rlo_ = new double[nrows0_] ;
275 CoinMemcpyN(rowLower,len,rlo_) ;
276
277 return ; }
278
279void CoinPrePostsolveMatrix::setRowUpper (const double *rowUpper, int lenParam)
280
281{ int len ;
282
283 if (lenParam < 0)
284 { len = nrows_ ; }
285 else
286 if (lenParam > nrows0_)
287 { throw CoinError("length exceeds allocated size",
288 "setRowUpper","CoinPrePostsolveMatrix") ; }
289 else
290 { len = lenParam ; }
291
292 if (rup_ == 0) rup_ = new double[nrows0_] ;
293 CoinMemcpyN(rowUpper,len,rup_) ;
294
295 return ; }
296
297void CoinPrePostsolveMatrix::setRowPrice (const double *rowSol, int lenParam)
298
299{ int len ;
300
301 if (lenParam < 0)
302 { len = nrows_ ; }
303 else
304 if (lenParam > nrows0_)
305 { throw CoinError("length exceeds allocated size",
306 "setRowPrice","CoinPrePostsolveMatrix") ; }
307 else
308 { len = lenParam ; }
309
310 if (rowduals_ == 0) rowduals_ = new double[nrows0_] ;
311 CoinMemcpyN(rowSol,len,rowduals_) ;
312
313 return ; }
314
315void CoinPrePostsolveMatrix::setRowActivity (const double *rowAct, int lenParam)
316
317{ int len ;
318
319 if (lenParam < 0)
320 { len = nrows_ ; }
321 else
322 if (lenParam > nrows0_)
323 { throw CoinError("length exceeds allocated size",
324 "setRowActivity","CoinPrePostsolveMatrix") ; }
325 else
326 { len = lenParam ; }
327
328 if (acts_ == 0) acts_ = new double[nrows0_] ;
329 CoinMemcpyN(rowAct,len,acts_) ;
330
331 return ; }
332
333
334
335/*
336 Methods to set the status vectors for a basis. Note that we need to allocate
337 colstat_ and rowstat_ as a single vector, to maintain compatibility with
338 OsiPresolve and ClpPresolve.
339
340 The `using ::getStatus' declaration is required to get the compiler to
341 consider the getStatus helper function defined in CoinWarmStartBasis.hpp.
342*/
343
344void CoinPrePostsolveMatrix::setStructuralStatus (const char *strucStatus,
345 int lenParam)
346
347{ int len ;
348 using ::getStatus ;
349
350 if (lenParam < 0)
351 { len = ncols_ ; }
352 else
353 if (lenParam > ncols0_)
354 { throw CoinError("length exceeds allocated size",
355 "setStructuralStatus","CoinPrePostsolveMatrix") ; }
356 else
357 { len = lenParam ; }
358
359 if (colstat_ == 0)
360 { colstat_ = new unsigned char[ncols0_+nrows0_] ;
361# ifdef ZEROFAULT
362 CoinZeroN(colstat_,ncols0_+nrows0_) ;
363# endif
364 rowstat_ = colstat_+ncols0_ ; }
365 for (int j = 0 ; j < len ; j++)
366 { Status statj = Status(getStatus(strucStatus,j)) ;
367 setColumnStatus(j,statj) ; }
368
369 return ; }
370
371
372void CoinPrePostsolveMatrix::setArtificialStatus (const char *artifStatus,
373 int lenParam)
374
375{ int len ;
376 using ::getStatus ;
377
378 if (lenParam < 0)
379 { len = nrows_ ; }
380 else
381 if (lenParam > nrows0_)
382 { throw CoinError("length exceeds allocated size",
383 "setArtificialStatus","CoinPrePostsolveMatrix") ; }
384 else
385 { len = lenParam ; }
386
387 if (colstat_ == 0)
388 { colstat_ = new unsigned char[ncols0_+nrows0_] ;
389# ifdef ZEROFAULT
390 CoinZeroN(colstat_,ncols0_+nrows0_) ;
391# endif
392 rowstat_ = colstat_+ncols0_ ; }
393 for (int i = 0 ; i < len ; i++)
394 { Status stati = Status(getStatus(artifStatus,i)) ;
395 setRowStatus(i,stati) ; }
396
397 return ; }
398
399/*
400 This routine initialises structural and artificial status given a
401 CoinWarmStartBasis as the parameter.
402*/
403
404void CoinPrePostsolveMatrix::setStatus (const CoinWarmStartBasis *basis)
405
406{ setStructuralStatus(basis->getStructuralStatus(),
407 basis->getNumStructural()) ;
408 setArtificialStatus(basis->getArtificialStatus(),
409 basis->getNumArtificial()) ;
410
411 return ; }
412
413/*
414 This routine returns structural and artificial status in the form of a
415 CoinWarmStartBasis object.
416
417 What to do when CoinPrePostsolveMatrix::Status == superBasic? There's
418 no analog in CoinWarmStartBasis::Status.
419*/
420
421CoinWarmStartBasis *CoinPrePostsolveMatrix::getStatus ()
422
423{ int n = ncols_ ;
424 int m = nrows_ ;
425 CoinWarmStartBasis *wsb = new CoinWarmStartBasis() ;
426 wsb->setSize(n,m) ;
427 for (int j = 0 ; j < n ; j++)
428 { CoinWarmStartBasis::Status statj =
429 CoinWarmStartBasis::Status(getColumnStatus(j)) ;
430 wsb->setStructStatus(j,statj) ; }
431 for (int i = 0 ; i < m ; i++)
432 { CoinWarmStartBasis::Status stati =
433 CoinWarmStartBasis::Status(getRowStatus(i)) ;
434 wsb->setArtifStatus(i,stati) ; }
435
436 return (wsb) ; }
437#endif
438/*
439 Set the status of a non-basic artificial variable based on the
440 variable's value and bounds.
441*/
442
443void CoinPrePostsolveMatrix::setRowStatusUsingValue (int iRow)
444
445{ double value = acts_[iRow];
446 double lower = rlo_[iRow];
447 double upper = rup_[iRow];
448 if (lower<-1.0e20&&upper>1.0e20) {
449 setRowStatus(iRow,isFree);
450 } else if (fabs(lower-value)<=ztolzb_) {
451 setRowStatus(iRow,atLowerBound);
452 } else if (fabs(upper-value)<=ztolzb_) {
453 setRowStatus(iRow,atUpperBound);
454 } else {
455 setRowStatus(iRow,superBasic);
456 }
457}
458
459/*
460 Set the status of a non-basic structural variable based on the
461 variable's value and bounds.
462*/
463
464void CoinPrePostsolveMatrix::setColumnStatusUsingValue(int iColumn)
465{
466 double value = sol_[iColumn];
467 double lower = clo_[iColumn];
468 double upper = cup_[iColumn];
469 if (lower<-1.0e20&&upper>1.0e20) {
470 setColumnStatus(iColumn,isFree);
471 } else if (fabs(lower-value)<=ztolzb_) {
472 setColumnStatus(iColumn,atLowerBound);
473 } else if (fabs(upper-value)<=ztolzb_) {
474 setColumnStatus(iColumn,atUpperBound);
475 } else {
476 setColumnStatus(iColumn,superBasic);
477 }
478}
479#ifndef SLIM_CLP
480
481
482/*
483 Simple routines to return a constant character string for the status value.
484 Separate row and column routines for convenience.
485*/
486
487const char *CoinPrePostsolveMatrix::columnStatusString (int j) const
488
489{ Status statj = getColumnStatus(j) ;
490
491 switch (statj)
492 { case isFree:
493 { return ("NBF") ; }
494 case basic:
495 { return ("B") ; }
496 case atUpperBound:
497 { return ("NBUB") ; }
498 case atLowerBound:
499 { return ("NBLB") ; }
500 case superBasic:
501 { return ("SB") ; }
502 default:
503 { return ("INVALID") ; } } }
504
505const char *CoinPrePostsolveMatrix::rowStatusString (int j) const
506
507{ Status statj = getRowStatus(j) ;
508
509 switch (statj)
510 { case isFree:
511 { return ("NBF") ; }
512 case basic:
513 { return ("B") ; }
514 case atUpperBound:
515 { return ("NBUB") ; }
516 case atLowerBound:
517 { return ("NBLB") ; }
518 case superBasic:
519 { return ("SB") ; }
520 default:
521 { return ("INVALID") ; } } }
522#endif
523