1/* $Id: CoinPresolveDoubleton.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#include <math.h>
8
9#include "CoinFinite.hpp"
10#include "CoinHelperFunctions.hpp"
11#include "CoinPresolveMatrix.hpp"
12
13#include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW
14#include "CoinPresolveZeros.hpp"
15#include "CoinPresolveFixed.hpp"
16#include "CoinPresolveDoubleton.hpp"
17
18#include "CoinPresolvePsdebug.hpp"
19#include "CoinMessage.hpp"
20
21#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
22#include "CoinPresolvePsdebug.hpp"
23#endif
24
25
26namespace { /* begin unnamed local namespace */
27
28/*
29 This routine does the grunt work needed to substitute x for y in all rows i
30 where coeff[i,y] != 0. We have
31
32 y = (c - a*x)/b = c/b + (-a/b)*x
33
34 Suppose we're fixing row i. We need to adjust the row bounds by
35 -coeff[i,y]*(c/b) and coeff[i,x] by coeff[i,y]*(-a/b). The value
36 c/b is passed as the bounds_factor, and -a/b as the coeff_factor.
37
38 row0 is the doubleton row. It is assumed that coeff[row0,y] has been
39 removed from the column major representation before this routine is
40 called. (Otherwise, we'd have to check for it to avoid a useless row
41 update.)
42
43 Both the row and col representations are updated. There are two cases:
44
45 * coeff[i,x] != 0:
46 in the column rep, modify coeff[i,x];
47 in the row rep, modify coeff[i,x] and drop coeff[i,y].
48
49 * coeff[i,x] == 0 (i.e., non-existent):
50 in the column rep, add coeff[i,x]; mcstrt is modified if the column
51 must be moved;
52 in the row rep, convert coeff[i,y] to coeff[i,x].
53
54 The row and column reps are inconsistent during the routine and at
55 completion. In the row rep, column x and y are updated except for
56 the doubleton row, and in the column rep only column x is updated
57 except for coeff[row0,x]. On return, column y and row row0 will be deleted
58 and consistency will be restored.
59*/
60
61 bool elim_doubleton (const char *
62#ifdef PRESOLVE_DEBUG
63msg
64#endif
65 ,
66 CoinBigIndex *mcstrt,
67 double *rlo, double *rup,
68 double *colels,
69 int *hrow, int *hcol,
70 int *hinrow, int *hincol,
71 presolvehlink *clink, int ncols,
72 CoinBigIndex *mrstrt, double *rowels,
73 double coeff_factor,
74 double bounds_factor,
75 int
76#ifdef PRESOLVE_DEBUG
77 row0
78#endif
79 , int icolx, int icoly)
80
81{
82 CoinBigIndex kcsx = mcstrt[icolx];
83 CoinBigIndex kcex = kcsx + hincol[icolx];
84
85# if PRESOLVE_DEBUG
86 printf("%s %d x=%d y=%d cf=%g bf=%g nx=%d yrows=(", msg,
87 row0, icolx, icoly, coeff_factor, bounds_factor, hincol[icolx]);
88# endif
89/*
90 Open a loop to scan column y. For each nonzero coefficient (row,y),
91 update column x and the row bounds for the row.
92
93 The initial assert checks that we're properly updating column x.
94*/
95 CoinBigIndex base = mcstrt[icoly];
96 int numberInY = hincol[icoly];
97 for (int kwhere = 0 ; kwhere < numberInY ; kwhere++)
98 { PRESOLVEASSERT(kcex == kcsx+hincol[icolx]) ;
99 CoinBigIndex kcoly = base+kwhere;
100/*
101 Look for coeff[row,x], then update accordingly.
102*/
103 double coeffy = colels[kcoly] ;
104 double delta = coeffy*coeff_factor ;
105 int row = hrow[kcoly] ;
106 CoinBigIndex kcolx = presolve_find_row1(row,kcsx,kcex,hrow) ;
107# if PRESOLVE_DEBUG
108 printf("%d%s ",row,(kcolx<kcex)?"+":"") ;
109# endif
110/*
111 Case 1: coeff[i,x] != 0: update it in column and row reps; drop coeff[i,y]
112 from row rep.
113*/
114 if (kcolx < kcex)
115 { colels[kcolx] += delta ;
116
117 CoinBigIndex kmi = presolve_find_col(icolx,mrstrt[row],
118 mrstrt[row]+hinrow[row],hcol) ;
119 rowels[kmi] = colels[kcolx] ;
120 presolve_delete_from_row(row,icoly,mrstrt,hinrow,hcol,rowels) ; }
121/*
122 Case 2: coeff[i,x] == 0: add it in the column rep; convert coeff[i,y] in
123 the row rep. presolve_expand_col ensures an empty entry exists at the
124 end of the column. The location of column x may change with expansion.
125*/
126 else
127 { bool no_mem = presolve_expand_col(mcstrt,colels,hrow,hincol,
128 clink,ncols,icolx) ;
129 if (no_mem)
130 return (true) ;
131
132 kcsx = mcstrt[icolx] ;
133 kcex = mcstrt[icolx]+hincol[icolx] ;
134 // recompute y as well
135 base = mcstrt[icoly];
136
137 hrow[kcex] = row ;
138 colels[kcex] = delta ;
139 hincol[icolx]++ ;
140 kcex++ ;
141
142 CoinBigIndex k2 = presolve_find_col(icoly,mrstrt[row],
143 mrstrt[row]+hinrow[row],hcol) ;
144 hcol[k2] = icolx ;
145 rowels[k2] = delta ; }
146/*
147 Update the row bounds, if necessary. Avoid updating finite infinity.
148*/
149 if (bounds_factor != 0.0)
150 { delta = coeffy*bounds_factor ;
151 if (-PRESOLVE_INF < rlo[row])
152 rlo[row] -= delta ;
153 if (rup[row] < PRESOLVE_INF)
154 rup[row] -= delta ; } }
155
156# if PRESOLVE_DEBUG
157 printf(")\n") ;
158# endif
159
160 return (false) ; }
161
162} /* end unnamed local namespace */
163
164
165/*
166 * It is always the case that one of the variables of a doubleton
167 * will be (implied) free, but neither will necessarily be a singleton.
168 * Since in the case of a doubleton the number of non-zero entries
169 * will never increase, though, it makes sense to always eliminate them.
170 *
171 * The col rep and row rep must be consistent.
172 */
173const CoinPresolveAction
174 *doubleton_action::presolve (CoinPresolveMatrix *prob,
175 const CoinPresolveAction *next)
176
177{
178 double startTime = 0.0;
179 int startEmptyRows=0;
180 int startEmptyColumns = 0;
181 if (prob->tuning_) {
182 startTime = CoinCpuTime();
183 startEmptyRows = prob->countEmptyRows();
184 startEmptyColumns = prob->countEmptyCols();
185 }
186 double *colels = prob->colels_;
187 int *hrow = prob->hrow_;
188 CoinBigIndex *mcstrt = prob->mcstrt_;
189 int *hincol = prob->hincol_;
190 int ncols = prob->ncols_;
191
192 double *clo = prob->clo_;
193 double *cup = prob->cup_;
194
195 double *rowels = prob->rowels_;
196 int *hcol = prob->hcol_;
197 CoinBigIndex *mrstrt = prob->mrstrt_;
198 int *hinrow = prob->hinrow_;
199 int nrows = prob->nrows_;
200
201 double *rlo = prob->rlo_;
202 double *rup = prob->rup_;
203
204 presolvehlink *clink = prob->clink_;
205 presolvehlink *rlink = prob->rlink_;
206
207 const unsigned char *integerType = prob->integerType_;
208
209 double *cost = prob->cost_;
210
211 int numberLook = prob->numberRowsToDo_;
212 int iLook;
213 int * look = prob->rowsToDo_;
214 const double ztolzb = prob->ztolzb_;
215
216 action * actions = new action [nrows];
217 int nactions = 0;
218
219 int *zeros = prob->usefulColumnInt_; //new int[ncols];
220 int nzeros = 0;
221
222 int *fixed = zeros+ncols; //new int[ncols];
223 int nfixed = 0;
224
225 unsigned char *rowstat = prob->rowstat_;
226 double *acts = prob->acts_;
227 double * sol = prob->sol_;
228
229 bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
230# if PRESOLVE_CONSISTENCY
231 presolve_consistent(prob) ;
232 presolve_links_ok(prob) ;
233# endif
234
235 // wasfor (int irow=0; irow<nrows; irow++)
236 for (iLook=0;iLook<numberLook;iLook++) {
237 int irow = look[iLook];
238 if (hinrow[irow] == 2 &&
239 fabs(rup[irow] - rlo[irow]) <= ZTOLDP) {
240 double rhs = rlo[irow];
241 CoinBigIndex krs = mrstrt[irow];
242 int icolx, icoly;
243 CoinBigIndex k;
244
245 icolx = hcol[krs];
246 icoly = hcol[krs+1];
247 if (hincol[icolx]<=0||hincol[icoly]<=0) {
248 // should never happen ?
249 //printf("JJF - doubleton column %d has %d entries and %d has %d\n",
250 // icolx,hincol[icolx],icoly,hincol[icoly]);
251 continue;
252 }
253 // check size
254 if (fabs(rowels[krs]) < ZTOLDP2 || fabs(rowels[krs+1]) < ZTOLDP2)
255 continue;
256 // See if prohibited for any reason
257 if (prob->colProhibited(icolx) || prob->colProhibited(icolx))
258 continue;
259
260 // don't bother with fixed variables
261 if (!(fabs(cup[icolx] - clo[icolx]) < ZTOLDP) &&
262 !(fabs(cup[icoly] - clo[icoly]) < ZTOLDP)) {
263 double coeffx, coeffy;
264 /* find this row in each of the columns */
265 CoinBigIndex krowx = presolve_find_row(irow, mcstrt[icolx], mcstrt[icolx] + hincol[icolx], hrow);
266 CoinBigIndex krowy = presolve_find_row(irow, mcstrt[icoly], mcstrt[icoly] + hincol[icoly], hrow);
267
268/*
269 Check for integrality: If one variable is integer, keep it and substitute
270 for the continuous variable. If both are integer, substitute only for the
271 forms x = k * y (k integral and non-empty intersection on bounds on x)
272 or x = 1-y, where both x and y are binary.
273
274 flag bits for integerStatus: 1>>0 x integer
275 1>>1 y integer
276*/
277 int integerStatus=0;
278 if (integerType[icolx]) {
279 if (integerType[icoly]) {
280 // both integer
281 int good = 0;
282 double rhs2 = rhs;
283 double value;
284 value=colels[krowx];
285 if (value<0.0) {
286 value = - value;
287 rhs2 += 1;
288 }
289 if (cup[icolx]==1.0&&clo[icolx]==0.0&&fabs(value-1.0)<1.0e-7)
290 good =1;
291 value=colels[krowy];
292 if (value<0.0) {
293 value = - value;
294 rhs2 += 1;
295 }
296 if (cup[icoly]==1.0&&clo[icoly]==0.0&&fabs(value-1.0)<1.0e-7)
297 good |= 2;
298 if (good==3&&fabs(rhs2-1.0)<1.0e-7)
299 integerStatus = 3;
300 else
301 integerStatus=-1;
302 if (integerStatus==-1&&!rhs) {
303 // maybe x = k * y;
304 double value1 = colels[krowx];
305 double value2 = colels[krowy];
306 double ratio;
307 bool swap=false;
308 if (fabs(value1)>fabs(value2)) {
309 ratio = value1/value2;
310 } else {
311 ratio = value2/value1;
312 swap=true;
313 }
314 ratio=fabs(ratio);
315 if (fabs(ratio-floor(ratio+0.5))<1.0e-12) {
316 // possible
317 integerStatus = swap ? 2 : 1;
318 //printf("poss type %d\n",integerStatus);
319 }
320 }
321 } else {
322 integerStatus = 1;
323 }
324 } else if (integerType[icoly]) {
325 integerStatus = 2;
326 }
327 if (integerStatus<0) {
328 // can still take in some cases
329 bool canDo=false;
330 double value1 = colels[krowx];
331 double value2 = colels[krowy];
332 double ratio;
333 bool swap=false;
334 double rhsRatio;
335 if (fabs(value1)>fabs(value2)) {
336 ratio = value1/value2;
337 rhsRatio = rhs/value1;
338 } else {
339 ratio = value2/value1;
340 rhsRatio = rhs/value2;
341 swap=true;
342 }
343 ratio=fabs(ratio);
344 if (fabs(ratio-floor(ratio+0.5))<1.0e-12) {
345 // possible
346 integerStatus = swap ? 2 : 1;
347 // but check rhs
348 if (rhsRatio==floor(rhsRatio+0.5))
349 canDo=true;
350 }
351#ifdef COIN_DEVELOP2
352 if (canDo)
353 printf("Good CoinPresolveDoubleton icolx %d (%g and bounds %g %g) icoly %d (%g and bound %g %g) - rhs %g\n",
354 icolx,colels[krowx],clo[icolx],cup[icolx],
355 icoly,colels[krowy],clo[icoly],cup[icoly],rhs);
356 else
357 printf("Bad CoinPresolveDoubleton icolx %d (%g) icoly %d (%g) - rhs %g\n",
358 icolx,colels[krowx],icoly,colels[krowy],rhs);
359#endif
360 if (!canDo)
361 continue;
362 }
363 if (integerStatus == 2) {
364 CoinSwap(icoly,icolx);
365 CoinSwap(krowy,krowx);
366 }
367
368 // HAVE TO JIB WITH ABOVE swapS
369 // if x's coefficient is something like 1000, but y's only something like -1,
370 // then when we postsolve, if x's is close to being out of tolerance,
371 // then y is very likely to be (because y==1000x) . (55)
372 // It it interesting that the number of doubletons found may depend
373 // on which column is substituted away (this is true of baxter.mps).
374 if (!integerStatus) {
375 if (fabs(colels[krowy]) < fabs(colels[krowx])) {
376 CoinSwap(icoly,icolx);
377 CoinSwap(krowy,krowx);
378 }
379 }
380
381#if 0
382 //?????
383 if (integerType[icolx] &&
384 clo[icoly] != -PRESOLVE_INF &&
385 cup[icoly] != PRESOLVE_INF) {
386 continue;
387 }
388#endif
389
390 {
391 CoinBigIndex kcs = mcstrt[icoly];
392 CoinBigIndex kce = kcs + hincol[icoly];
393 for (k=kcs; k<kce; k++) {
394 if (hinrow[hrow[k]] == 1) {
395 break;
396 }
397 }
398 // let singleton rows be taken care of first
399 if (k<kce)
400 continue;
401 }
402
403 coeffx = colels[krowx];
404 coeffy = colels[krowy];
405
406 // it is possible that both x and y are singleton columns
407 // that can cause problems
408 if (hincol[icolx] == 1 && hincol[icoly] == 1)
409 continue;
410
411 // BE CAUTIOUS and avoid very large relative differences
412 // if this is not done in baxter, then the computed solution isn't optimal,
413 // but gets it in 11995 iterations; the postsolve goes to iteration 16181.
414 // with this, the solution is optimal, but takes 18825 iters; postsolve 18871.
415#if 0
416 if (fabs(coeffx) * max_coeff_factor <= fabs(coeffy))
417 continue;
418#endif
419
420#if 0
421 if (only_zero_rhs && rhs != 0)
422 continue;
423
424 if (reject_doubleton(mcstrt, colels, hrow, hincol,
425 -coeffx / coeffy,
426 max_coeff_ratio,
427 irow, icolx, icoly))
428 continue;
429#endif
430
431 // common equations are of the form ax + by = 0, or x + y >= lo
432 {
433 PRESOLVE_DETAIL_PRINT(printf("pre_doubleton %dC %dC %dR E\n",
434 icoly,icolx,irow));
435 action *s = &actions[nactions];
436 nactions++;
437
438 s->row = irow;
439 s->icolx = icolx;
440
441 s->clox = clo[icolx];
442 s->cupx = cup[icolx];
443 s->costx = cost[icolx];
444
445 s->icoly = icoly;
446 s->costy = cost[icoly];
447
448 s->rlo = rlo[irow];
449
450 s->coeffx = coeffx;
451
452 s->coeffy = coeffy;
453
454 s->ncolx = hincol[icolx];
455
456 s->ncoly = hincol[icoly];
457 if (s->ncoly<s->ncolx) {
458 // Take out row
459 s->colel = presolve_dupmajor(colels,hrow,hincol[icoly],
460 mcstrt[icoly],irow) ;
461 s->ncolx=0;
462 } else {
463 s->colel = presolve_dupmajor(colels,hrow,hincol[icolx],
464 mcstrt[icolx],irow) ;
465 s->ncoly=0;
466 }
467 }
468
469 /*
470 * This moves the bounds information for y onto x,
471 * making y free and allowing us to substitute it away.
472 *
473 * a x + b y = c
474 * l1 <= x <= u1
475 * l2 <= y <= u2 ==>
476 *
477 * l2 <= (c - a x) / b <= u2
478 * b/-a > 0 ==> (b l2 - c) / -a <= x <= (b u2 - c) / -a
479 * b/-a < 0 ==> (b u2 - c) / -a <= x <= (b l2 - c) / -a
480 */
481 {
482 double lo1 = -PRESOLVE_INF;
483 double up1 = PRESOLVE_INF;
484
485 //PRESOLVEASSERT((coeffx < 0) == (coeffy/-coeffx < 0));
486 // (coeffy/-coeffx < 0) == (coeffy<0 == coeffx<0)
487 if (-PRESOLVE_INF < clo[icoly]) {
488 if (coeffx * coeffy < 0)
489 lo1 = (coeffy * clo[icoly] - rhs) / -coeffx;
490 else
491 up1 = (coeffy * clo[icoly] - rhs) / -coeffx;
492 }
493
494 if (cup[icoly] < PRESOLVE_INF) {
495 if (coeffx * coeffy < 0)
496 up1 = (coeffy * cup[icoly] - rhs) / -coeffx;
497 else
498 lo1 = (coeffy * cup[icoly] - rhs) / -coeffx;
499 }
500
501 // costy y = costy ((c - a x) / b) = (costy c)/b + x (costy -a)/b
502 // the effect of maxmin cancels out
503 cost[icolx] += cost[icoly] * (-coeffx / coeffy);
504
505 prob->change_bias(cost[icoly] * rhs / coeffy);
506
507 if (0 /*integerType[icolx]*/) {
508 abort();
509 /* no change possible for now */
510#if 0
511 lo1 = trunc(lo1);
512 up1 = trunc(up1);
513
514 /* trunc(3.5) == 3.0 */
515 /* trunc(-3.5) == -3.0 */
516
517 /* I think this is ok */
518 if (lo1 > clo[icolx]) {
519 (clo[icolx] <= 0.0)
520 clo[icolx] = ? ilo
521
522 clo[icolx] = ilo;
523 cup[icolx] = iup;
524 }
525#endif
526 } else {
527 double lo2 = CoinMax(clo[icolx], lo1);
528 double up2 = CoinMin(cup[icolx], up1);
529 if (lo2 > up2) {
530 if (lo2 <= up2 + prob->feasibilityTolerance_||fixInfeasibility) {
531 // If close to integer then go there
532 double nearest = floor(lo2+0.5);
533 if (fabs(nearest-lo2)<2.0*prob->feasibilityTolerance_) {
534 lo2 = nearest;
535 up2 = nearest;
536 } else {
537 lo2 = up2;
538 }
539 } else {
540 prob->status_ |= 1;
541 prob->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
542 prob->messages())
543 <<icolx
544 <<lo2
545 <<up2
546 <<CoinMessageEol;
547 break;
548 }
549 }
550 clo[icolx] = lo2;
551 cup[icolx] = up2;
552
553 if (rowstat&&sol) {
554 // update solution and basis
555 int basisChoice=0;
556 int numberBasic=0;
557 double movement = 0 ;
558 if (prob->columnIsBasic(icolx))
559 numberBasic++;
560 if (prob->columnIsBasic(icoly))
561 numberBasic++;
562 if (prob->rowIsBasic(irow))
563 numberBasic++;
564 if (sol[icolx]<=lo2+ztolzb) {
565 movement = lo2-sol[icolx] ;
566 sol[icolx] = lo2;
567 prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::atLowerBound);
568 } else if (sol[icolx]>=up2-ztolzb) {
569 movement = up2-sol[icolx] ;
570 sol[icolx] = up2;
571 prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::atUpperBound);
572 } else {
573 basisChoice=1;
574 }
575 if (numberBasic>1)
576 prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::basic);
577/*
578 We need to compensate if x was forced to move. Beyond that, even if x didn't
579 move, we've forced y = (c-ax)/b, and that might not have been true before. So
580 even if x didn't move, y may have moved. Note that the constant term c/b is
581 subtracted out as the constraints are modified, so we don't include it when
582 calculating movement for y.
583*/
584 if (movement)
585 { CoinBigIndex k;
586 for (k = mcstrt[icolx] ; k < mcstrt[icolx]+hincol[icolx] ; k++)
587 { int row = hrow[k];
588 if (hinrow[row])
589 acts[row] += movement*colels[k]; } }
590 movement = (-coeffx*sol[icolx]/coeffy)-sol[icoly] ;
591 if (movement)
592 { for (k = mcstrt[icoly] ;
593 k < mcstrt[icoly]+hincol[icoly] ;
594 k++)
595 { int row = hrow[k];
596 if (hinrow[row])
597 acts[row] += movement*colels[k]; } }
598 }
599 if (lo2 == up2)
600 fixed[nfixed++] = icolx;
601 }
602 }
603
604 // Update next set of actions
605 {
606 prob->addCol(icolx);
607 int i,kcs,kce;
608 kcs = mcstrt[icoly];
609 kce = kcs + hincol[icoly];
610 for (i=kcs;i<kce;i++) {
611 int row = hrow[i];
612 prob->addRow(row);
613 }
614 kcs = mcstrt[icolx];
615 kce = kcs + hincol[icolx];
616 for (i=kcs;i<kce;i++) {
617 int row = hrow[i];
618 prob->addRow(row);
619 }
620 }
621
622/*
623 Empty irow in the column-major matrix. Deleting the coefficient for
624 (irow,icoly) is a bit costly (given that we're about to drop the whole
625 column), but saves the trouble of checking for it in elim_doubleton.
626*/
627 presolve_delete_from_col(irow,icolx,mcstrt,hincol,hrow,colels) ;
628 presolve_delete_from_col(irow,icoly,mcstrt,hincol,hrow,colels) ;
629/*
630 Drop irow in the row-major representation: set the length to 0
631 and reclaim the major vector space in bulk storage.
632*/
633 hinrow[irow] = 0;
634 PRESOLVE_REMOVE_LINK(rlink,irow);
635
636 /* transfer the colx factors to coly */
637 bool no_mem = elim_doubleton("ELIMD",
638 mcstrt, rlo, rup, colels,
639 hrow, hcol, hinrow, hincol,
640 clink, ncols,
641 mrstrt, rowels,
642 -coeffx / coeffy,
643 rhs / coeffy,
644 irow, icolx, icoly);
645 if (no_mem)
646 throwCoinError("out of memory",
647 "doubleton_action::presolve");
648
649
650 // eliminate coly entirely from the col rep
651 hincol[icoly] = 0;
652 PRESOLVE_REMOVE_LINK(clink, icoly);
653 cost[icoly] = 0.0;
654
655 rlo[irow] = 0.0;
656 rup[irow] = 0.0;
657
658 zeros[nzeros++] = icolx; // check for zeros
659
660 // strictly speaking, since we didn't adjust {clo,cup}[icoly]
661 // or {rlo,rup}[irow], this col/row may be infeasible,
662 // because the solution/activity would be 0, whereas the
663 // bounds may be non-zero.
664 }
665
666# if PRESOLVE_CONSISTENCY
667 presolve_consistent(prob) ;
668 presolve_links_ok(prob) ;
669# endif
670 }
671 }
672
673 if (nactions) {
674# if PRESOLVE_SUMMARY
675 printf("NDOUBLETONS: %d\n", nactions);
676# endif
677 action *actions1 = new action[nactions];
678 CoinMemcpyN(actions, nactions, actions1);
679
680 next = new doubleton_action(nactions, actions1, next);
681
682 if (nzeros)
683 next = drop_zero_coefficients_action::presolve(prob, zeros, nzeros, next);
684 if (nfixed)
685 next = remove_fixed_action::presolve(prob, fixed, nfixed, next);
686 }
687
688 //delete[]zeros;
689 //delete[]fixed;
690 deleteAction(actions,action*);
691
692 if (prob->tuning_) {
693 double thisTime=CoinCpuTime();
694 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
695 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
696 printf("CoinPresolveDoubleton(4) - %d rows, %d columns dropped in time %g, total %g\n",
697 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
698 }
699 return (next);
700}
701
702/*
703 Reintroduce the column (y) and doubleton row (irow) removed in presolve.
704 Correct the other column (x) involved in the doubleton, update the solution,
705 etc.
706
707 A fair amount of complication arises because the presolve transform saves the
708 shorter of x or y. Postsolve thus includes portions to restore either.
709*/
710void doubleton_action::postsolve(CoinPostsolveMatrix *prob) const
711{
712 const action *const actions = actions_;
713 const int nactions = nactions_;
714
715 double *colels = prob->colels_;
716 int *hrow = prob->hrow_;
717 CoinBigIndex *mcstrt = prob->mcstrt_;
718 int *hincol = prob->hincol_;
719 int *link = prob->link_;
720
721 double *clo = prob->clo_;
722 double *cup = prob->cup_;
723
724 double *rlo = prob->rlo_;
725 double *rup = prob->rup_;
726
727 double *dcost = prob->cost_;
728
729 double *sol = prob->sol_;
730 double *acts = prob->acts_;
731 double *rowduals = prob->rowduals_;
732 double *rcosts = prob->rcosts_;
733
734 unsigned char *colstat = prob->colstat_;
735 unsigned char *rowstat = prob->rowstat_;
736
737 const double maxmin = prob->maxmin_;
738
739# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
740 char *cdone = prob->cdone_;
741 char *rdone = prob->rdone_;
742# endif
743
744 CoinBigIndex &free_list = prob->free_list_;
745
746 const double ztolzb = prob->ztolzb_;
747 const double ztoldj = prob->ztoldj_;
748
749 int nrows = prob->nrows_;
750 // Arrays to rebuild the unsaved column.
751 int * index1 = new int[nrows];
752 double * element1 = new double[nrows];
753 CoinZeroN(element1,nrows);
754
755# if PRESOLVE_CONSISTENCY
756 presolve_check_threads(prob) ;
757# endif
758# if PRESOLVE_DEBUG
759 presolve_check_sol(prob) ;
760# endif
761/*
762 The outer loop: step through the doubletons in this array of actions.
763 The first activity is to unpack the doubleton.
764*/
765 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
766
767 int irow = f->row;
768 double lo0 = f->clox;
769 double up0 = f->cupx;
770
771
772 double coeffx = f->coeffx;
773 double coeffy = f->coeffy;
774 int jcolx = f->icolx;
775 int jcoly = f->icoly;
776
777 double rhs = f->rlo;
778/*
779 jcolx is in the problem (for whatever reason), and the doubleton row (irow)
780 and column (jcoly) have only been processed by empty row/column postsolve
781 (i.e., reintroduced with length 0).
782*/
783 PRESOLVEASSERT(cdone[jcolx] && rdone[irow]==DROP_ROW);
784 PRESOLVEASSERT(cdone[jcoly]==DROP_COL);
785
786/*
787 Restore bounds for doubleton row, bounds and objective coefficient for x,
788 objective for y.
789
790 Original comment: restoration of rlo and rup likely isn't necessary.
791*/
792 rlo[irow] = f->rlo;
793 rup[irow] = f->rlo;
794
795 clo[jcolx] = lo0;
796 cup[jcolx] = up0;
797
798 dcost[jcolx] = f->costx;
799 dcost[jcoly] = f->costy;
800
801#if PRESOLVE_DEBUG
802/* Original comment: I've forgotten what this is about
803
804 Loss of significant digits through cancellation, with possible inflation
805 when divided by coeffy below? -- lh, 040831 --
806*/
807 if ((rhs < 0) == ((coeffx * sol[jcolx]) < 0) &&
808 fabs(rhs - coeffx * sol[jcolx]) * 100 < rhs &&
809 fabs(rhs - coeffx * sol[jcolx]) * 100 < (coeffx * sol[jcolx]))
810 printf("DANGEROUS RHS??? %g %g %g\n",
811 rhs, coeffx * sol[jcolx],
812 (rhs - coeffx * sol[jcolx]));
813#endif
814/*
815 Set primal solution for y (including status) and row activity for the
816 doubleton row. The motivation (up in presolve) for wanting coeffx < coeffy
817 is to avoid inflation into sol[y]. Since this is a (satisfied) equality,
818 activity is the rhs value and the logical is nonbasic.
819*/
820 sol[jcoly] = (rhs-coeffx*sol[jcolx])/coeffy;
821 acts[irow] = rhs;
822 if (rowstat)
823 prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound);
824/*
825 Time to get into the correction/restoration of coefficients for columns x
826 and y, with attendant correction of row bounds and activities. Accumulate
827 partial reduced costs (missing the contribution from the doubleton row) so
828 that we can eventually calculate a dual for the doubleton row.
829*/
830 double djy = maxmin * dcost[jcoly];
831 double djx = maxmin * dcost[jcolx];
832 double bounds_factor = rhs/coeffy;
833/*
834 We saved column y in the action, so we'll use it to reconstruct column x.
835 There are two aspects: correction of existing x coefficients, and fill in.
836 Given
837 coeffx'[k] = coeffx[k]+coeffy[k]*coeff_factor
838 we have
839 coeffx[k] = coeffx'[k]-coeffy[k]*coeff_factor
840 where
841 coeff_factor = -coeffx[dblton]/coeffy[dblton].
842
843 Keep in mind that the major vector stored in the action does not include
844 the coefficient from the doubleton row --- the doubleton coefficients are
845 held in coeffx and coeffy.
846*/
847 if (f->ncoly) {
848 int ncoly=f->ncoly-1; // as row taken out
849 double multiplier = coeffx/coeffy;
850 //printf("Current colx %d\n",jcolx);
851 int * indy = reinterpret_cast<int *>(f->colel+ncoly);
852/*
853 Rebuild a threaded column y, starting with the end of the thread and working
854 back to the beginning. In the process, accumulate corrections to column x
855 in element1 and index1. Fix row bounds and activity as we go (add back the
856 constant correction removed in presolve), and accumulate contributions to
857 the reduced cost for y.
858
859 The PRESOLVEASSERT says this row should already be present.
860*/
861 int ystart = NO_LINK;
862 int nX=0;
863 int i,iRow;
864 for (i=0; i<ncoly; ++i) {
865 int iRow = indy[i];
866 PRESOLVEASSERT(rdone[iRow]);
867
868 double yValue = f->colel[i];
869
870 // undo elim_doubleton(1)
871 if (-PRESOLVE_INF < rlo[iRow])
872 rlo[iRow] += yValue * bounds_factor;
873
874 // undo elim_doubleton(2)
875 if (rup[iRow] < PRESOLVE_INF)
876 rup[iRow] += yValue * bounds_factor;
877
878 acts[iRow] += yValue * bounds_factor;
879
880 djy -= rowduals[iRow] * yValue;
881/*
882 Link the coefficient into column y: Acquire the first free slot in the
883 bulk arrays and store the row index and coefficient. Then link the slot
884 in front of coefficients we've already processed.
885*/
886 CoinBigIndex k = free_list;
887 assert(k >= 0 && k < prob->bulk0_) ;
888 free_list = link[free_list];
889 hrow[k] = iRow;
890 colels[k] = yValue;
891 link[k] = ystart;
892 ystart = k;
893/*
894 Calculate and store the correction to the x coefficient.
895*/
896 yValue *= multiplier;
897 element1[iRow]=yValue;
898 index1[nX++]=iRow;
899 }
900# if PRESOLVE_CONSISTENCY
901 presolve_check_free_list(prob) ;
902# endif
903/*
904 Handle the coefficients of the doubleton row.
905*/
906 {
907 double yValue = coeffy;
908
909 CoinBigIndex k = free_list;
910 assert(k >= 0 && k < prob->bulk0_) ;
911 free_list = link[free_list];
912 hrow[k] = irow;
913 colels[k] = yValue;
914 link[k] = ystart;
915 ystart = k;
916
917 yValue *= multiplier;
918 element1[irow]=yValue;
919 index1[nX++]=irow;
920 }
921/*
922 Attach the threaded column y to mcstrt and record the length.
923*/
924 mcstrt[jcoly] = ystart;
925 hincol[jcoly] = f->ncoly;
926/*
927 Now integrate the corrections to column x. First thing to do is find the
928 end of the column. While we're doing that, correct any existing entries.
929 This complicates life because the correction could cancel the existing
930 coefficient and we don't want to leave an explicit zero. In this case we
931 relink the column around it. (last is a little misleading --- it's actually
932 `last nonzero'. If we haven't seen a nonzero yet, the relink goes to mcstrt.)
933 The freed slot is linked at the beginning of the free list.
934*/
935 CoinBigIndex k=mcstrt[jcolx];
936 CoinBigIndex last = NO_LINK;
937 int numberInColumn = hincol[jcolx];
938 int numberToDo=numberInColumn;
939 for (i=0; i<numberToDo; ++i) {
940 iRow = hrow[k];
941 assert (iRow>=0&&iRow<nrows);
942 double value = colels[k]+element1[iRow];
943 element1[iRow]=0.0;
944 if (fabs(value)>=1.0e-15) {
945 colels[k]=value;
946 last=k;
947 k = link[k];
948 if (iRow != irow)
949 djx -= rowduals[iRow] * value;
950 } else {
951 numberInColumn--;
952 // add to free list
953 int nextk = link[k];
954 assert(free_list>=0);
955 link[k]=free_list;
956 free_list=k;
957 assert (k>=0);
958 k=nextk;
959 if (last!=NO_LINK)
960 link[last]=k;
961 else
962 mcstrt[jcolx]=k;
963 }
964 }
965/*
966 We've found the end of column x. Any remaining nonzeros in element1 will be
967 fill in, which we link at the end of the column thread.
968*/
969 for (i=0;i<nX;i++) {
970 int iRow = index1[i];
971 double xValue = element1[iRow];
972 element1[iRow]=0.0;
973 if (fabs(xValue)>=1.0e-15) {
974 if (iRow != irow)
975 djx -= rowduals[iRow] * xValue;
976 numberInColumn++;
977 CoinBigIndex k = free_list;
978 assert(k >= 0 && k < prob->bulk0_) ;
979 free_list = link[free_list];
980 hrow[k] = iRow;
981 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
982 colels[k] = xValue;
983 if (last!=NO_LINK)
984 link[last] = k;
985 else
986 mcstrt[jcolx]=k;
987 last = k;
988 }
989 }
990
991# if PRESOLVE_CONSISTENCY
992 presolve_check_free_list(prob) ;
993# endif
994
995/*
996 Whew! Tidy up column x and we're done.
997*/
998 link[last]=NO_LINK;
999 assert(numberInColumn);
1000 hincol[jcolx] = numberInColumn;
1001/*
1002 Of course, we could have saved column x in the action. Now we need to
1003 regenerate coefficients of column y.
1004 Given
1005 coeffx'[k] = coeffx[k]+coeffy[k]*coeff_factor
1006 we have
1007 coeffy[k] = (coeffx'[k]-coeffx[k])*(1/coeff_factor)
1008 where
1009 coeff_factor = -coeffx[dblton]/coeffy[dblton].
1010*/
1011 } else {
1012 int ncolx=f->ncolx-1;
1013 double multiplier = -coeffy/coeffx;
1014 int * indx = reinterpret_cast<int *> (f->colel+ncolx);
1015 //printf("Current colx %d\n",jcolx);
1016/*
1017 Scan existing column x to find the end. While we're at it, accumulate part
1018 of the new y coefficients in index1 and element1.
1019*/
1020 CoinBigIndex k=mcstrt[jcolx];
1021 int nX=0;
1022 int i,iRow;
1023 for (i=0; i<hincol[jcolx]-1; ++i) {
1024 if (colels[k]) {
1025 iRow = hrow[k];
1026 index1[nX++]=iRow;
1027 element1[iRow]=multiplier*colels[k];
1028 }
1029 k = link[k];
1030 }
1031 if (colels[k]) {
1032 iRow = hrow[k];
1033 index1[nX++]=iRow;
1034 element1[iRow]=multiplier*colels[k];
1035 }
1036/*
1037 Replace column x with the the original column x held in the doubleton
1038 action. We first move column x to the free list, then thread a column with
1039 the original coefficients, back to front. While we're at it, add the
1040 second part of the y coefficients to index1 and element1.
1041*/
1042 multiplier = - multiplier;
1043 link[k] = free_list;
1044 free_list = mcstrt[jcolx];
1045 int xstart = NO_LINK;
1046 for (i=0; i<ncolx; ++i) {
1047 int iRow = indx[i];
1048 PRESOLVEASSERT(rdone[iRow]);
1049
1050 double xValue = f->colel[i];
1051 //printf("x %d %d %g\n",i,indx[i],f->colel[i]);
1052 CoinBigIndex k = free_list;
1053 assert(k >= 0 && k < prob->bulk0_) ;
1054 free_list = link[free_list];
1055 hrow[k] = iRow;
1056 colels[k] = xValue;
1057 link[k] = xstart;
1058 xstart = k;
1059
1060 djx -= rowduals[iRow] * xValue;
1061
1062 xValue *= multiplier;
1063 if (!element1[iRow]) {
1064 element1[iRow]=xValue;
1065 index1[nX++]=iRow;
1066 } else {
1067 element1[iRow]+=xValue;
1068 }
1069 }
1070# if PRESOLVE_CONSISTENCY
1071 presolve_check_free_list(prob) ;
1072# endif
1073/*
1074 The same, for the doubleton row.
1075*/
1076 {
1077 double xValue = coeffx;
1078 CoinBigIndex k = free_list;
1079 assert(k >= 0 && k < prob->bulk0_) ;
1080 free_list = link[free_list];
1081 hrow[k] = irow;
1082 colels[k] = xValue;
1083 link[k] = xstart;
1084 xstart = k;
1085
1086 xValue *= multiplier;
1087 if (!element1[irow]) {
1088 element1[irow]=xValue;
1089 index1[nX++]=irow;
1090 } else {
1091 element1[irow]+=xValue;
1092 }
1093 }
1094/*
1095 Link the new column x to mcstrt and set the length.
1096*/
1097 mcstrt[jcolx] = xstart;
1098 hincol[jcolx] = f->ncolx;
1099/*
1100 Now get to work building a threaded column y from the nonzeros in element1.
1101 As before, build the thread in reverse.
1102*/
1103 int ystart = NO_LINK;
1104 int n=0;
1105 for (i=0;i<nX;i++) {
1106 int iRow = index1[i];
1107 PRESOLVEASSERT(rdone[iRow] || iRow == irow);
1108 double yValue = element1[iRow];
1109 element1[iRow]=0.0;
1110 if (fabs(yValue)>=1.0e-12) {
1111 n++;
1112 CoinBigIndex k = free_list;
1113 assert(k >= 0 && k < prob->bulk0_) ;
1114 free_list = link[free_list];
1115 hrow[k] = iRow;
1116 colels[k] = yValue;
1117 link[k] = ystart;
1118 ystart = k;
1119 }
1120 }
1121# if PRESOLVE_CONSISTENCY
1122 presolve_check_free_list(prob) ;
1123# endif
1124/*
1125 Tidy up --- link the new column into mcstrt and set the length.
1126*/
1127 mcstrt[jcoly] = ystart;
1128 assert(n);
1129 hincol[jcoly] = n;
1130/*
1131 Now that we have the original y, we can scan it and do the corrections to
1132 the row bounds and activity, and get a start on a reduced cost for y.
1133*/
1134 k = mcstrt[jcoly];
1135 int ny = hincol[jcoly];
1136 for (i=0; i<ny; ++i) {
1137 int row = hrow[k];
1138 double coeff = colels[k];
1139 k = link[k];
1140
1141 if (row != irow) {
1142
1143 // undo elim_doubleton(1)
1144 if (-PRESOLVE_INF < rlo[row])
1145 rlo[row] += coeff * bounds_factor;
1146
1147 // undo elim_doubleton(2)
1148 if (rup[row] < PRESOLVE_INF)
1149 rup[row] += coeff * bounds_factor;
1150
1151 acts[row] += coeff * bounds_factor;
1152
1153 djy -= rowduals[row] * coeff;
1154 }
1155 }
1156/*
1157 Scan the new column x and calculate reduced costs. This could be integrated
1158 into the previous section where the original column x is restored.
1159
1160 ok --- let's try it, then.
1161
1162 k = mcstrt[jcolx];
1163 int nx = hincol[jcolx];
1164
1165 for ( i=0; i<nx; ++i) {
1166 int row = hrow[k];
1167 double coeff = colels[k];
1168 k = link[k];
1169
1170 if (row != irow) {
1171 djx -= rowduals[row] * coeff;
1172 }
1173 }
1174*/
1175 }
1176/*
1177 Sanity? The only assignment to coeffx is f->coeffx! Ditto for coeffy.
1178*/
1179 assert (fabs(coeffx-f->coeffx)<1.0e-6&&fabs(coeffy-f->coeffy)<1.0e-6);
1180/*
1181 Time to calculate a dual for the doubleton row, and settle the status of x
1182 and y. Ideally, we'll leave x at whatever nonbasic status it currently has
1183 and make y basic. There's a potential problem, however: Remember that we
1184 transferred bounds from y to x when we eliminated y. If those bounds were
1185 tighter than x's original bounds, we may not be able to maintain x at its
1186 present status, or even as nonbasic.
1187
1188 We'll make two claims here:
1189
1190 * If the dual value for the doubleton row is chosen to keep the reduced
1191 cost djx of col x at its prior value, then the reduced cost djy of col
1192 y will be 0. (Crank through the linear algebra to convince yourself.)
1193
1194 * If the bounds on x have loosened, then it must be possible to make y
1195 nonbasic, because we've transferred the tight bound back to y. (Yeah,
1196 I'm waving my hands. But it sounds good. -- lh, 040907 --)
1197
1198 So ... if we can maintain x nonbasic, then we need to set y basic, which
1199 means we should calculate rowduals[dblton] so that rcost[jcoly] == 0. We
1200 may need to change the status of x (an artifact of loosening a bound when
1201 x was previously a fixed variable).
1202
1203 If we need to push x into the basis, then we calculate rowduals[dblton] so
1204 that rcost[jcolx] == 0 and make y nonbasic.
1205*/
1206 // printf("djs x - %g (%g), y - %g (%g)\n",djx,coeffx,djy,coeffy);
1207 if (colstat)
1208 { bool basicx = prob->columnIsBasic(jcolx) ;
1209 bool nblbxok = (fabs(lo0 - sol[jcolx]) < ztolzb) &&
1210 (rcosts[jcolx] >= -ztoldj) ;
1211 bool nbubxok = (fabs(up0 - sol[jcolx]) < ztolzb) &&
1212 (rcosts[jcolx] <= ztoldj) ;
1213 if (basicx || nblbxok || nbubxok)
1214 { if (!basicx)
1215 { if (nblbxok)
1216 { prob->setColumnStatus(jcolx,
1217 CoinPrePostsolveMatrix::atLowerBound) ; }
1218 else
1219 if (nbubxok)
1220 { prob->setColumnStatus(jcolx,
1221 CoinPrePostsolveMatrix::atUpperBound) ; } }
1222 prob->setColumnStatus(jcoly,CoinPrePostsolveMatrix::basic);
1223 rowduals[irow] = djy / coeffy;
1224 rcosts[jcolx] = djx - rowduals[irow] * coeffx;
1225#if 0
1226 if (prob->columnIsBasic(jcolx))
1227 assert (fabs(rcosts[jcolx])<1.0e-5);
1228#endif
1229 rcosts[jcoly] = 0.0;
1230 } else {
1231 prob->setColumnStatus(jcolx,CoinPrePostsolveMatrix::basic);
1232 prob->setColumnStatusUsingValue(jcoly);
1233 rowduals[irow] = djx / coeffx;
1234 rcosts[jcoly] = djy - rowduals[irow] * coeffy;
1235 rcosts[jcolx] = 0.0;
1236 }
1237 } else {
1238 // No status array
1239 // this is the coefficient we need to force col y's reduced cost to 0.0;
1240 // for example, this is obviously true if y is a singleton column
1241 rowduals[irow] = djy / coeffy;
1242 rcosts[jcoly] = 0.0;
1243 }
1244
1245# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1246/*
1247 Mark the column and row as processed by doubleton action. The check integrity
1248 of the threaded matrix.
1249*/
1250 cdone[jcoly] = DOUBLETON;
1251 rdone[irow] = DOUBLETON;
1252 presolve_check_threads(prob) ;
1253#endif
1254# if PRESOLVE_DEBUG
1255/*
1256 Confirm accuracy of reduced cost for columns x and y.
1257*/
1258 {
1259 CoinBigIndex k = mcstrt[jcolx];
1260 int nx = hincol[jcolx];
1261 double dj = maxmin * dcost[jcolx];
1262
1263 for (int i=0; i<nx; ++i) {
1264 int row = hrow[k];
1265 double coeff = colels[k];
1266 k = link[k];
1267
1268 dj -= rowduals[row] * coeff;
1269 }
1270 if (! (fabs(rcosts[jcolx] - dj) < 100*ZTOLDP))
1271 printf("BAD DOUBLE X DJ: %d %d %g %g\n",
1272 irow, jcolx, rcosts[jcolx], dj);
1273 rcosts[jcolx]=dj;
1274 }
1275 {
1276 CoinBigIndex k = mcstrt[jcoly];
1277 int ny = hincol[jcoly];
1278 double dj = maxmin * dcost[jcoly];
1279
1280 for (int i=0; i<ny; ++i) {
1281 int row = hrow[k];
1282 double coeff = colels[k];
1283 k = link[k];
1284
1285 dj -= rowduals[row] * coeff;
1286 //printf("b %d coeff %g dual %g dj %g\n",
1287 // row,coeff,rowduals[row],dj);
1288 }
1289 if (! (fabs(rcosts[jcoly] - dj) < 100*ZTOLDP))
1290 printf("BAD DOUBLE Y DJ: %d %d %g %g\n",
1291 irow, jcoly, rcosts[jcoly], dj);
1292 rcosts[jcoly]=dj;
1293 }
1294# endif
1295 }
1296/*
1297 Done at last. Delete the scratch arrays.
1298*/
1299
1300 delete [] index1;
1301 delete [] element1;
1302}
1303
1304
1305doubleton_action::~doubleton_action()
1306{
1307 for (int i=nactions_-1; i>=0; i--) {
1308 delete[]actions_[i].colel;
1309 }
1310 deleteAction(actions_,action*);
1311}
1312
1313
1314
1315static double *doubleton_mult;
1316static int *doubleton_id;
1317void check_doubletons(const CoinPresolveAction * paction)
1318{
1319 const CoinPresolveAction * paction0 = paction;
1320
1321 if (paction) {
1322 check_doubletons(paction->next);
1323
1324 if (strcmp(paction0->name(), "doubleton_action") == 0) {
1325 const doubleton_action *daction =
1326 reinterpret_cast<const doubleton_action *>(paction0);
1327 for (int i=daction->nactions_-1; i>=0; --i) {
1328 int icolx = daction->actions_[i].icolx;
1329 int icoly = daction->actions_[i].icoly;
1330 double coeffx = daction->actions_[i].coeffx;
1331 double coeffy = daction->actions_[i].coeffy;
1332
1333 doubleton_mult[icoly] = -coeffx/coeffy;
1334 doubleton_id[icoly] = icolx;
1335 }
1336 }
1337 }
1338}
1339
1340#if PRESOLVE_DEBUG
1341void check_doubletons1(const CoinPresolveAction * paction,
1342 int ncols)
1343#else
1344void check_doubletons1(const CoinPresolveAction * /*paction*/,
1345 int /*ncols*/)
1346#endif
1347{
1348#if PRESOLVE_DEBUG
1349 doubleton_mult = new double[ncols];
1350 doubleton_id = new int[ncols];
1351 int i;
1352 for ( i=0; i<ncols; ++i)
1353 doubleton_id[i] = i;
1354 check_doubletons(paction);
1355 double minmult = 1.0;
1356 int minid = -1;
1357 for ( i=0; i<ncols; ++i) {
1358 double mult = 1.0;
1359 int j = i;
1360 if (doubleton_id[j] != j) {
1361 printf("MULTS (%d): ", j);
1362 while (doubleton_id[j] != j) {
1363 printf("%d %g, ", doubleton_id[j], doubleton_mult[j]);
1364 mult *= doubleton_mult[j];
1365 j = doubleton_id[j];
1366 }
1367 printf(" == %g\n", mult);
1368 if (minmult > fabs(mult)) {
1369 minmult = fabs(mult);
1370 minid = i;
1371 }
1372 }
1373 }
1374 if (minid != -1)
1375 printf("MIN MULT: %d %g\n", minid, minmult);
1376#endif
1377}
1378