1/* $Id: IdiSolve.cpp 1753 2011-06-19 16:27:26Z 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 "CoinPragma.hpp"
7#include <stdio.h>
8#include <stdarg.h>
9#include <stdlib.h>
10#include <math.h>
11#include "CoinHelperFunctions.hpp"
12#include "Idiot.hpp"
13#define FIT
14#ifdef FIT
15#define HISTORY 8
16#else
17#define HISTORY 7
18#endif
19#define NSOLVE HISTORY-1
20static void solveSmall(int nsolve, double **aIn, double **a, double * b)
21{
22 int i, j;
23 /* copy */
24 for (i = 0; i < nsolve; i++) {
25 for (j = 0; j < nsolve; j++) {
26 a[i][j] = aIn[i][j];
27 }
28 }
29 for (i = 0; i < nsolve; i++) {
30 /* update using all previous */
31 double diagonal;
32 int j;
33 for (j = i; j < nsolve; j++) {
34 int k;
35 double value = a[i][j];
36 for (k = 0; k < i; k++) {
37 value -= a[k][i] * a[k][j];
38 }
39 a[i][j] = value;
40 }
41 diagonal = a[i][i];
42 if (diagonal < 1.0e-20) {
43 diagonal = 0.0;
44 } else {
45 diagonal = 1.0 / sqrt(diagonal);
46 }
47 a[i][i] = diagonal;
48 for (j = i + 1; j < nsolve; j++) {
49 a[i][j] *= diagonal;
50 }
51 }
52 /* update */
53 for (i = 0; i < nsolve; i++) {
54 int j;
55 double value = b[i];
56 for (j = 0; j < i; j++) {
57 value -= b[j] * a[j][i];
58 }
59 value *= a[i][i];
60 b[i] = value;
61 }
62 for (i = nsolve - 1; i >= 0; i--) {
63 int j;
64 double value = b[i];
65 for (j = i + 1; j < nsolve; j++) {
66 value -= b[j] * a[i][j];
67 }
68 value *= a[i][i];
69 b[i] = value;
70 }
71}
72IdiotResult
73Idiot::objval(int nrows, int ncols, double * rowsol , double * colsol,
74 double * pi, double * /*djs*/, const double * cost ,
75 const double * /*rowlower*/,
76 const double * rowupper, const double * /*lower*/,
77 const double * /*upper*/, const double * elemnt,
78 const int * row, const CoinBigIndex * columnStart,
79 const int * length, int extraBlock, int * rowExtra,
80 double * solExtra, double * elemExtra, double * /*upperExtra*/,
81 double * costExtra, double weight)
82{
83 IdiotResult result;
84 double objvalue = 0.0;
85 double sum1 = 0.0, sum2 = 0.0;
86 int i;
87 for (i = 0; i < nrows; i++) {
88 rowsol[i] = -rowupper[i];
89 }
90 for (i = 0; i < ncols; i++) {
91 CoinBigIndex j;
92 double value = colsol[i];
93 if (value) {
94 objvalue += value * cost[i];
95 if (elemnt) {
96 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
97 int irow = row[j];
98 rowsol[irow] += elemnt[j] * value;
99 }
100 } else {
101 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
102 int irow = row[j];
103 rowsol[irow] += value;
104 }
105 }
106 }
107 }
108 /* adjust to make as feasible as possible */
109 /* no */
110 if (extraBlock) {
111 for (i = 0; i < extraBlock; i++) {
112 double element = elemExtra[i];
113 int irow = rowExtra[i];
114 objvalue += solExtra[i] * costExtra[i];
115 rowsol[irow] += solExtra[i] * element;
116 }
117 }
118 for (i = 0; i < nrows; i++) {
119 double value = rowsol[i];
120 sum1 += fabs(value);
121 sum2 += value * value;
122 pi[i] = -2.0 * weight * value;
123 }
124 result.infeas = sum1;
125 result.objval = objvalue;
126 result.weighted = objvalue + weight * sum2;
127 result.sumSquared = sum2;
128 return result;
129}
130IdiotResult
131Idiot::IdiSolve(
132 int nrows, int ncols, double * rowsol , double * colsol,
133 double * pi, double * djs, const double * origcost , double * rowlower,
134 double * rowupper, const double * lower,
135 const double * upper, const double * elemnt,
136 const int * row, const CoinBigIndex * columnStart,
137 const int * length, double * lambda,
138 int maxIts, double mu, double drop,
139 double maxmin, double offset,
140 int strategy, double djTol, double djExit, double djFlag,
141 CoinThreadRandom * randomNumberGenerator)
142{
143 IdiotResult result;
144 int i, j, k, iter;
145 double value = 0.0, objvalue = 0.0, weightedObj = 0.0;
146 double tolerance = 1.0e-8;
147 double * history[HISTORY+1];
148 int ncolx;
149 int nChange;
150 int extraBlock = 0;
151 int * rowExtra = NULL;
152 double * solExtra = NULL;
153 double * elemExtra = NULL;
154 double * upperExtra = NULL;
155 double * costExtra = NULL;
156 double * useCostExtra = NULL;
157 double * saveExtra = NULL;
158 double * cost = NULL;
159 double saveValue = 1.0e30;
160 double saveOffset = offset;
161 double useOffset = offset;
162 /*#define NULLVECTOR*/
163#ifndef NULLVECTOR
164 int nsolve = NSOLVE;
165#else
166 int nsolve = NSOLVE + 1; /* allow for null vector */
167#endif
168 int nflagged;
169 double * thetaX;
170 double * djX;
171 double * bX;
172 double * vX;
173 double ** aX;
174 double **aworkX;
175 double ** allsum;
176 double * saveSol = 0;
177 const double * useCost = cost;
178 double bestSol = 1.0e60;
179 double weight = 0.5 / mu;
180 char * statusSave = new char[2*ncols];
181 char * statusWork = statusSave + ncols;
182#define DJTEST 5
183 double djSave[DJTEST];
184 double largestDj = 0.0;
185 double smallestDj = 1.0e60;
186 double maxDj = 0.0;
187 int doFull = 0;
188#define SAVEHISTORY 10
189#define EVERY (2*SAVEHISTORY)
190#define AFTER SAVEHISTORY*(HISTORY+1)
191#define DROP 5
192 double after = AFTER;
193 double obj[DROP];
194 double kbad = 0, kgood = 0;
195 if (strategy & 128) after = 999999; /* no acceleration at all */
196 for (i = 0; i < DROP; i++) {
197 obj[i] = 1.0e70;
198 }
199 //#define TWO_GOES
200#ifdef TWO_GOES
201 double * pi2 = new double [nrows];
202 double * rowsol2 = new double [nrows];
203#endif
204 allsum = new double * [nsolve];
205 aX = new double * [nsolve];
206 aworkX = new double * [nsolve];
207 thetaX = new double[nsolve];
208 vX = new double[nsolve];
209 bX = new double[nsolve];
210 djX = new double[nsolve];
211 allsum[0] = pi;
212 for (i = 0; i < nsolve; i++) {
213 if (i) allsum[i] = new double[nrows];
214 aX[i] = new double[nsolve];
215 aworkX[i] = new double[nsolve];
216 }
217 /* check = rows */
218 for (i = 0; i < nrows; i++) {
219 if (rowupper[i] - rowlower[i] > tolerance) {
220 extraBlock++;
221 }
222 }
223 cost = new double[ncols];
224 memset(rowsol, 0, nrows * sizeof(double));
225 for (i = 0; i < ncols; i++) {
226 CoinBigIndex j;
227 double value = origcost[i] * maxmin;
228 double value2 = colsol[i];
229 if (elemnt) {
230 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
231 int irow = row[j];
232 value += elemnt[j] * lambda[irow];
233 rowsol[irow] += elemnt[j] * value2;
234 }
235 } else {
236 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
237 int irow = row[j];
238 value += lambda[irow];
239 rowsol[irow] += value2;
240 }
241 }
242 cost[i] = value;
243 }
244 if (extraBlock) {
245 rowExtra = new int[extraBlock];
246 solExtra = new double[extraBlock];
247 elemExtra = new double[extraBlock];
248 upperExtra = new double[extraBlock];
249 costExtra = new double[extraBlock];
250 saveExtra = new double[extraBlock];
251 extraBlock = 0;
252 int nbad = 0;
253 for (i = 0; i < nrows; i++) {
254 if (rowupper[i] - rowlower[i] > tolerance) {
255 double smaller, difference;
256 double value;
257 saveExtra[extraBlock] = rowupper[i];
258 if (fabs(rowupper[i]) > fabs(rowlower[i])) {
259 smaller = rowlower[i];
260 value = -1.0;
261 } else {
262 smaller = rowupper[i];
263 value = 1.0;
264 }
265 if (fabs(smaller) > 1.0e10) {
266 if (!nbad)
267 COIN_DETAIL_PRINT(printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
268 i, smaller));
269 nbad++;
270 if (rowupper[i] < 0.0 || rowlower[i] > 0.0)
271 abort();
272 if (fabs(rowupper[i]) > fabs(rowlower[i])) {
273 rowlower[i] = -0.9e10;
274 smaller = rowlower[i];
275 value = -1.0;
276 } else {
277 rowupper[i] = 0.9e10;
278 saveExtra[extraBlock] = rowupper[i];
279 smaller = rowupper[i];
280 value = 1.0;
281 }
282 }
283 difference = rowupper[i] - rowlower[i];
284 difference = CoinMin(difference, 1.0e31);
285 rowupper[i] = smaller;
286 elemExtra[extraBlock] = value;
287 solExtra[extraBlock] = (rowupper[i] - rowsol[i]) / value;
288 if (solExtra[extraBlock] < 0.0) solExtra[extraBlock] = 0.0;
289 if (solExtra[extraBlock] > difference) solExtra[extraBlock] = difference;
290 costExtra[extraBlock] = lambda[i] * value;
291 upperExtra[extraBlock] = difference;
292 rowsol[i] += value * solExtra[extraBlock];
293 rowExtra[extraBlock++] = i;
294 }
295 }
296 if (nbad)
297 COIN_DETAIL_PRINT(printf("%d bad values - results may be wrong\n", nbad));
298 }
299 for (i = 0; i < nrows; i++) {
300 offset += lambda[i] * rowsol[i];
301 }
302 if ((strategy & 256) != 0) {
303 /* save best solution */
304 saveSol = new double[ncols];
305 CoinMemcpyN(colsol, ncols, saveSol);
306 if (extraBlock) {
307 useCostExtra = new double[extraBlock];
308 memset(useCostExtra, 0, extraBlock * sizeof(double));
309 }
310 useCost = origcost;
311 useOffset = saveOffset;
312 } else {
313 useCostExtra = costExtra;
314 useCost = cost;
315 useOffset = offset;
316 }
317 ncolx = ncols + extraBlock;
318 for (i = 0; i < HISTORY + 1; i++) {
319 history[i] = new double[ncolx];
320 }
321 for (i = 0; i < DJTEST; i++) {
322 djSave[i] = 1.0e30;
323 }
324 for (i = 0; i < ncols; i++) {
325 if (upper[i] - lower[i]) {
326 statusSave[i] = 0;
327 } else {
328 statusSave[i] = 1;
329 }
330 }
331 // for two pass method
332 int start[2];
333 int stop[2];
334 int direction = -1;
335 start[0] = 0;
336 stop[0] = ncols;
337 start[1] = 0;
338 stop[1] = 0;
339 iter = 0;
340 for (; iter < maxIts; iter++) {
341 double sum1 = 0.0, sum2 = 0.0, djtot;
342 double lastObj = 1.0e70;
343 int good = 0, doScale = 0;
344 if (strategy & 16) {
345 int ii = iter / EVERY + 1;
346 ii = ii * EVERY;
347 if (iter > ii - HISTORY * 2 && (iter & 1) == 0) {
348 double * x = history[HISTORY-1];
349 for (i = HISTORY - 1; i > 0; i--) {
350 history[i] = history[i-1];
351 }
352 history[0] = x;
353 CoinMemcpyN(colsol, ncols, history[0]);
354 CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
355 }
356 }
357 if ((iter % SAVEHISTORY) == 0 || doFull) {
358 if ((strategy & 16) == 0) {
359 double * x = history[HISTORY-1];
360 for (i = HISTORY - 1; i > 0; i--) {
361 history[i] = history[i-1];
362 }
363 history[0] = x;
364 CoinMemcpyN(colsol, ncols, history[0]);
365 CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
366 }
367 }
368 /* start full try */
369 if ((iter % EVERY) == 0 || doFull) {
370 // for next pass
371 direction = - direction;
372 // randomize.
373 // The cast is to avoid gcc compiler warning
374 int kcol = static_cast<int>(ncols * randomNumberGenerator->randomDouble());
375 if (kcol == ncols)
376 kcol = ncols - 1;
377 if (direction > 0) {
378 start[0] = kcol;
379 stop[0] = ncols;
380 start[1] = 0;
381 stop[1] = kcol;
382 } else {
383 start[0] = kcol;
384 stop[0] = -1;
385 start[1] = ncols - 1;
386 stop[1] = kcol;
387 }
388 int itry = 0;
389 /*if ((strategy&16)==0) {
390 double * x=history[HISTORY-1];
391 for (i=HISTORY-1;i>0;i--) {
392 history[i]=history[i-1];
393 }
394 history[0]=x;
395 CoinMemcpyN(colsol,ncols,history[0]);
396 CoinMemcpyN(solExtra,extraBlock,history[0]+ncols);
397 }*/
398 while (!good) {
399 itry++;
400#define MAXTRY 5
401 if (iter > after && doScale < 2 && itry < MAXTRY) {
402 /* now full one */
403 for (i = 0; i < nrows; i++) {
404 rowsol[i] = -rowupper[i];
405 }
406 sum2 = 0.0;
407 objvalue = 0.0;
408 memset(pi, 0, nrows * sizeof(double));
409 djtot = 0.0;
410 {
411 double * theta = thetaX;
412 double * dj = djX;
413 double * b = bX;
414 double ** a = aX;
415 double ** awork = aworkX;
416 double * v = vX;
417 double c;
418#ifdef FIT
419 int ntot = 0, nsign = 0, ngood = 0, mgood[4] = {0, 0, 0, 0};
420 double diff1, diff2, val0, val1, val2, newValue;
421 CoinMemcpyN(colsol, ncols, history[HISTORY-1]);
422 CoinMemcpyN(solExtra, extraBlock, history[HISTORY-1] + ncols);
423#endif
424 dj[0] = 0.0;
425 for (i = 1; i < nsolve; i++) {
426 dj[i] = 0.0;
427 memset(allsum[i], 0, nrows * sizeof(double));
428 }
429 for (i = 0; i < ncols; i++) {
430 double value2 = colsol[i];
431 if (value2 > lower[i] + tolerance) {
432 if(value2 < (upper[i] - tolerance)) {
433 int k;
434 objvalue += value2 * cost[i];
435#ifdef FIT
436 ntot++;
437 val0 = history[0][i];
438 val1 = history[1][i];
439 val2 = history[2][i];
440 diff1 = val0 - val1;
441 diff2 = val1 - val2;
442 if (diff1*diff2 >= 0.0) {
443 nsign++;
444 if (fabs(diff1) < fabs(diff2)) {
445 int ii = static_cast<int>(fabs(4.0 * diff1 / diff2));
446 if (ii == 4) ii = 3;
447 mgood[ii]++;
448 ngood++;
449 }
450 if (fabs(diff1) < 0.75 * fabs(diff2)) {
451 newValue = val1 + (diff1 * diff2) / (diff2 - diff1);
452 } else {
453 newValue = val1 + 4.0 * diff1;
454 }
455 } else {
456 newValue = 0.333333333 * (val0 + val1 + val2);
457 }
458 if (newValue > upper[i] - tolerance) {
459 newValue = upper[i];
460 } else if (newValue < lower[i] + tolerance) {
461 newValue = lower[i];
462 }
463 history[HISTORY-1][i] = newValue;
464#endif
465 for (k = 0; k < HISTORY - 1; k++) {
466 value = history[k][i] - history[k+1][i];
467 dj[k] += value * cost[i];
468 v[k] = value;
469 }
470 if (elemnt) {
471 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
472 int irow = row[j];
473 for (k = 0; k < HISTORY - 1; k++) {
474 allsum[k][irow] += elemnt[j] * v[k];
475 }
476 rowsol[irow] += elemnt[j] * value2;
477 }
478 } else {
479 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
480 int irow = row[j];
481 for (k = 0; k < HISTORY - 1; k++) {
482 allsum[k][irow] += v[k];
483 }
484 rowsol[irow] += value2;
485 }
486 }
487 } else {
488 /* at ub */
489 colsol[i] = upper[i];
490 value2 = colsol[i];
491 objvalue += value2 * cost[i];
492 if (elemnt) {
493 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
494 int irow = row[j];
495 rowsol[irow] += elemnt[j] * value2;
496 }
497 } else {
498 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
499 int irow = row[j];
500 rowsol[irow] += value2;
501 }
502 }
503 }
504 } else {
505 /* at lb */
506 if (value2) {
507 objvalue += value2 * cost[i];
508 if (elemnt) {
509 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
510 int irow = row[j];
511 rowsol[irow] += elemnt[j] * value2;
512 }
513 } else {
514 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
515 int irow = row[j];
516 rowsol[irow] += value2;
517 }
518 }
519 }
520 }
521 }
522#ifdef FIT
523 /*printf("total %d, same sign %d, good %d %d %d %d %d\n",
524 ntot,nsign,ngood,mgood[0],mgood[1],mgood[2],mgood[3]);*/
525#endif
526 if (extraBlock) {
527 for (i = 0; i < extraBlock; i++) {
528 double element = elemExtra[i];
529 int irow = rowExtra[i];
530 objvalue += solExtra[i] * costExtra[i];
531 if (solExtra[i] > tolerance
532 && solExtra[i] < (upperExtra[i] - tolerance)) {
533 double value2 = solExtra[i];
534 int k;
535 for (k = 0; k < HISTORY - 1; k++) {
536 value = history[k][i+ncols] - history[k+1][i+ncols];
537 dj[k] += value * costExtra[i];
538 allsum[k][irow] += element * value;
539 }
540 rowsol[irow] += element * value2;
541 } else {
542 double value2 = solExtra[i];
543 double element = elemExtra[i];
544 int irow = rowExtra[i];
545 rowsol[irow] += element * value2;
546 }
547 }
548 }
549#ifdef NULLVECTOR
550 if ((strategy & 64)) {
551 double djVal = dj[0];
552 for (i = 0; i < ncols - nrows; i++) {
553 double value2 = colsol[i];
554 if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
555 value = history[0][i] - history[1][i];
556 } else {
557 value = 0.0;
558 }
559 history[HISTORY][i] = value;
560 }
561 for (; i < ncols; i++) {
562 double value2 = colsol[i];
563 double delta;
564 int irow = i - (ncols - nrows);
565 double oldSum = allsum[0][irow];
566 if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
567 delta = history[0][i] - history[1][i];
568 } else {
569 delta = 0.0;
570 }
571 djVal -= delta * cost[i];
572 oldSum -= delta;
573 delta = - oldSum;
574 djVal += delta * cost[i];
575 history[HISTORY][i] = delta;
576 }
577 dj[HISTORY-1] = djVal;
578 djVal = 0.0;
579 for (i = 0; i < ncols; i++) {
580 double value2 = colsol[i];
581 if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance ||
582 i >= ncols - nrows) {
583 int k;
584 value = history[HISTORY][i];
585 djVal += value * cost[i];
586 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
587 int irow = row[j];
588 allsum[nsolve-1][irow] += value;
589 }
590 }
591 }
592 printf("djs %g %g\n", dj[HISTORY-1], djVal);
593 }
594#endif
595 for (i = 0; i < nsolve; i++) {
596 int j;
597 b[i] = 0.0;
598 for (j = 0; j < nsolve; j++) {
599 a[i][j] = 0.0;
600 }
601 }
602 c = 0.0;
603 for (i = 0; i < nrows; i++) {
604 double value = rowsol[i];
605 for (k = 0; k < nsolve; k++) {
606 v[k] = allsum[k][i];
607 b[k] += v[k] * value;
608 }
609 c += value * value;
610 for (k = 0; k < nsolve; k++) {
611 for (j = k; j < nsolve; j++) {
612 a[k][j] += v[k] * v[j];
613 }
614 }
615 }
616 sum2 = c;
617 if (itry == 1) {
618 lastObj = objvalue + weight * sum2;
619 }
620 for (k = 0; k < nsolve; k++) {
621 b[k] = - (weight * b[k] + 0.5 * dj[k]);
622 for (j = k; j < nsolve; j++) {
623 a[k][j] *= weight;
624 a[j][k] = a[k][j];
625 }
626 }
627 c *= weight;
628 for (k = 0; k < nsolve; k++) {
629 theta[k] = b[k];
630 }
631 solveSmall(nsolve, a, awork, theta);
632 if ((strategy & 64) != 0) {
633 value = 10.0;
634 for (k = 0; k < nsolve; k++) {
635 value = CoinMax(value, fabs(theta[k]));
636 }
637 if (value > 10.0 && ((logLevel_ & 4) != 0)) {
638 printf("theta %g %g %g\n", theta[0], theta[1], theta[2]);
639 }
640 value = 10.0 / value;
641 for (k = 0; k < nsolve; k++) {
642 theta[k] *= value;
643 }
644 }
645 for (i = 0; i < ncolx; i++) {
646 double valueh = 0.0;
647 for (k = 0; k < HISTORY - 1; k++) {
648 value = history[k][i] - history[k+1][i];
649 valueh += value * theta[k];
650 }
651#ifdef NULLVECTOR
652 value = history[HISTORY][i];
653 valueh += value * theta[HISTORY-1];
654#endif
655 history[HISTORY][i] = valueh;
656 }
657 }
658#ifdef NULLVECTOR
659 if ((strategy & 64)) {
660 for (i = 0; i < ncols - nrows; i++) {
661 if (colsol[i] <= lower[i] + tolerance
662 || colsol[i] >= (upper[i] - tolerance)) {
663 history[HISTORY][i] = 0.0;
664 }
665 }
666 tolerance = -tolerance; /* switch off test */
667 }
668#endif
669 if (!doScale) {
670 for (i = 0; i < ncols; i++) {
671 if (colsol[i] > lower[i] + tolerance
672 && colsol[i] < (upper[i] - tolerance)) {
673 value = history[HISTORY][i];
674 colsol[i] += value;
675 if (colsol[i] < lower[i] + tolerance) {
676 colsol[i] = lower[i];
677 } else if (colsol[i] > upper[i] - tolerance) {
678 colsol[i] = upper[i];
679 }
680 }
681 }
682 if (extraBlock) {
683 for (i = 0; i < extraBlock; i++) {
684 if (solExtra[i] > tolerance
685 && solExtra[i] < (upperExtra[i] - tolerance)) {
686 value = history[HISTORY][i+ncols];
687 solExtra[i] += value;
688 if (solExtra[i] < 0.0) {
689 solExtra[i] = 0.0;
690 } else if (solExtra[i] > upperExtra[i]) {
691 solExtra[i] = upperExtra[i];
692 }
693 }
694 }
695 }
696 } else {
697 double theta = 1.0;
698 double saveTheta = theta;
699 for (i = 0; i < ncols; i++) {
700 if (colsol[i] > lower[i] + tolerance
701 && colsol[i] < (upper[i] - tolerance)) {
702 value = history[HISTORY][i];
703 if (value > 0) {
704 if (theta * value + colsol[i] > upper[i]) {
705 theta = (upper[i] - colsol[i]) / value;
706 }
707 } else if (value < 0) {
708 if (colsol[i] + theta * value < lower[i]) {
709 theta = (lower[i] - colsol[i]) / value;
710 }
711 }
712 }
713 }
714 if (extraBlock) {
715 for (i = 0; i < extraBlock; i++) {
716 if (solExtra[i] > tolerance
717 && solExtra[i] < (upperExtra[i] - tolerance)) {
718 value = history[HISTORY][i+ncols];
719 if (value > 0) {
720 if (theta * value + solExtra[i] > upperExtra[i]) {
721 theta = (upperExtra[i] - solExtra[i]) / value;
722 }
723 } else if (value < 0) {
724 if (solExtra[i] + theta * value < 0.0) {
725 theta = -solExtra[i] / value;
726 }
727 }
728 }
729 }
730 }
731 if ((iter % 100 == 0) && (logLevel_ & 8) != 0) {
732 if (theta < saveTheta) {
733 printf(" - modified theta %g\n", theta);
734 }
735 }
736 for (i = 0; i < ncols; i++) {
737 if (colsol[i] > lower[i] + tolerance
738 && colsol[i] < (upper[i] - tolerance)) {
739 value = history[HISTORY][i];
740 colsol[i] += value * theta;
741 }
742 }
743 if (extraBlock) {
744 for (i = 0; i < extraBlock; i++) {
745 if (solExtra[i] > tolerance
746 && solExtra[i] < (upperExtra[i] - tolerance)) {
747 value = history[HISTORY][i+ncols];
748 solExtra[i] += value * theta;
749 }
750 }
751 }
752 }
753#ifdef NULLVECTOR
754 tolerance = fabs(tolerance); /* switch back on */
755#endif
756 if ((iter % 100) == 0 && (logLevel_ & 8) != 0) {
757 printf("\n");
758 }
759 }
760 good = 1;
761 result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
762 rowlower, rowupper, lower, upper,
763 elemnt, row, columnStart, length, extraBlock, rowExtra,
764 solExtra, elemExtra, upperExtra, useCostExtra,
765 weight);
766 weightedObj = result.weighted;
767 if (!iter) saveValue = weightedObj;
768 objvalue = result.objval;
769 sum1 = result.infeas;
770 if (saveSol) {
771 if (result.weighted < bestSol) {
772 COIN_DETAIL_PRINT(printf("%d %g better than %g\n", iter,
773 result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
774 bestSol = result.weighted;
775 CoinMemcpyN(colsol, ncols, saveSol);
776 }
777 }
778#ifdef FITz
779 if (iter > after) {
780 IdiotResult result2;
781 double ww, oo, ss;
782 if (extraBlock) abort();
783 result2 = objval(nrows, ncols, row2, sol2, pi2, djs, cost,
784 rowlower, rowupper, lower, upper,
785 elemnt, row, columnStart, extraBlock, rowExtra,
786 solExtra, elemExtra, upperExtra, costExtra,
787 weight);
788 ww = result2.weighted;
789 oo = result2.objval;
790 ss = result2.infeas;
791 printf("wobj %g obj %g inf %g last %g\n", ww, oo, ss, lastObj);
792 if (ww < weightedObj && ww < lastObj) {
793 printf(" taken");
794 ntaken++;
795 saving += weightedObj - ww;
796 weightedObj = ww;
797 objvalue = oo;
798 sum1 = ss;
799 CoinMemcpyN(row2, nrows, rowsol);
800 CoinMemcpyN(pi2, nrows, pi);
801 CoinMemcpyN(sol2, ncols, colsol);
802 result = objval(nrows, ncols, rowsol, colsol, pi, djs, cost,
803 rowlower, rowupper, lower, upper,
804 elemnt, row, columnStart, extraBlock, rowExtra,
805 solExtra, elemExtra, upperExtra, costExtra,
806 weight);
807 weightedObj = result.weighted;
808 objvalue = result.objval;
809 sum1 = result.infeas;
810 if (ww < weightedObj) abort();
811 } else {
812 printf(" not taken");
813 nottaken++;
814 }
815 }
816#endif
817 /*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/
818 if (weightedObj > lastObj + 1.0e-4 && itry < MAXTRY) {
819 if((logLevel_ & 16) != 0 && doScale) {
820 printf("Weighted objective from %g to %g **** bad move\n",
821 lastObj, weightedObj);
822 }
823 if (doScale) {
824 good = 1;
825 }
826 if ((strategy & 3) == 1) {
827 good = 0;
828 if (weightedObj > lastObj + djExit) {
829 if ((logLevel_ & 16) != 0) {
830 printf("Weighted objective from %g to %g ?\n", lastObj, weightedObj);
831 }
832 CoinMemcpyN(history[0], ncols, colsol);
833 CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
834 good = 1;
835 }
836 } else if ((strategy & 3) == 2) {
837 if (weightedObj > lastObj + 0.1 * maxDj) {
838 CoinMemcpyN(history[0], ncols, colsol);
839 CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
840 doScale++;
841 good = 0;
842 }
843 } else if ((strategy & 3) == 3) {
844 if (weightedObj > lastObj + 0.001 * maxDj) {
845 /*doScale++;*/
846 good = 0;
847 }
848 }
849 }
850 }
851 if ((iter % checkFrequency_) == 0) {
852 double best = weightedObj;
853 double test = obj[0];
854 for (i = 1; i < DROP; i++) {
855 obj[i-1] = obj[i];
856 if (best > obj[i]) best = obj[i];
857 }
858 obj[DROP-1] = best;
859 if (test - best < drop && (strategy & 8) == 0) {
860 if ((logLevel_ & 8) != 0) {
861 printf("Exiting as drop in %d its is %g after %d iterations\n",
862 DROP * checkFrequency_, test - best, iter);
863 }
864 goto RETURN;
865 }
866 }
867 if ((iter % logFreq_) == 0) {
868 double piSum = 0.0;
869 for (i = 0; i < nrows; i++) {
870 piSum += (rowsol[i] + rowupper[i]) * pi[i];
871 }
872 if ((logLevel_ & 2) != 0) {
873 printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n",
874 iter, sum1, objvalue * maxmin - useOffset, weightedObj - useOffset,
875 piSum * maxmin - useOffset, maxDj);
876 }
877 }
878 CoinMemcpyN(statusSave, ncols, statusWork);
879 nflagged = 0;
880 }
881 nChange = 0;
882 doFull = 0;
883 maxDj = 0.0;
884 // go through forwards or backwards and starting at odd places
885 int itry;
886#ifdef TWO_GOES
887 memcpy(pi2, pi, nrows * sizeof(double));
888 memcpy(rowsol2, rowsol, nrows * sizeof(double));
889#endif
890 for (itry = 0; itry < 2; itry++) {
891 int icol = start[itry];
892 int istop = stop[itry];
893#ifdef TWO_GOES
894 for (int iPar = 0; iPar < 2; iPar++) {
895 double * temp = pi;
896 pi = pi2;
897 pi2 = temp;
898 temp = rowsol;
899 rowsol = rowsol2;
900 rowsol2 = temp;
901 if (iPar == 0) {
902 istop = (icol + istop) >> 1;
903 } else {
904 icol = istop;
905 istop = stop[itry];
906 }
907#endif
908 for (; icol != istop; icol += direction) {
909 if (!statusWork[icol]) {
910 CoinBigIndex j;
911 double value = colsol[icol];
912 double djval = cost[icol];
913 double djval2, value2;
914 double theta, a, b, c;
915 if (elemnt) {
916 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
917 int irow = row[j];
918 djval -= elemnt[j] * pi[irow];
919 }
920 } else {
921 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
922 int irow = row[j];
923 djval -= pi[irow];
924 }
925 }
926 /*printf("xx iter %d seq %d djval %g value %g\n",
927 iter,i,djval,value);*/
928 if (djval > 1.0e-5) {
929 value2 = (lower[icol] - value);
930 } else {
931 value2 = (upper[icol] - value);
932 }
933 djval2 = djval * value2;
934 djval = fabs(djval);
935 if (djval > djTol) {
936 if (djval2 < -1.0e-4) {
937 double valuep, thetap;
938 nChange++;
939 if (djval > maxDj) maxDj = djval;
940 /*if (djval>3.55e6) {
941 printf("big\n");
942 }*/
943 a = 0.0;
944 b = 0.0;
945 c = 0.0;
946 djval2 = cost[icol];
947 if (elemnt) {
948 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
949 int irow = row[j];
950 double value = rowsol[irow];
951 c += value * value;
952 a += elemnt[j] * elemnt[j];
953 b += value * elemnt[j];
954 }
955 } else {
956 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
957 int irow = row[j];
958 double value = rowsol[irow];
959 c += value * value;
960 a += 1.0;
961 b += value;
962 }
963 }
964 a *= weight;
965 b = b * weight + 0.5 * djval2;
966 c *= weight;
967 /* solve */
968 theta = -b / a;
969 if ((strategy & 4) != 0) {
970 value2 = a * theta * theta + 2.0 * b * theta;
971 thetap = 2.0 * theta;
972 valuep = a * thetap * thetap + 2.0 * b * thetap;
973 if (valuep < value2 + djTol) {
974 theta = thetap;
975 kgood++;
976 } else {
977 kbad++;
978 }
979 }
980 if (theta > 0.0) {
981 if (theta < upper[icol] - colsol[icol]) {
982 value2 = theta;
983 } else {
984 value2 = upper[icol] - colsol[icol];
985 }
986 } else {
987 if (theta > lower[icol] - colsol[icol]) {
988 value2 = theta;
989 } else {
990 value2 = lower[icol] - colsol[icol];
991 }
992 }
993 colsol[icol] += value2;
994 objvalue += cost[icol] * value2;
995 if (elemnt) {
996 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
997 int irow = row[j];
998 double value;
999 rowsol[irow] += elemnt[j] * value2;
1000 value = rowsol[irow];
1001 pi[irow] = -2.0 * weight * value;
1002 }
1003 } else {
1004 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
1005 int irow = row[j];
1006 double value;
1007 rowsol[irow] += value2;
1008 value = rowsol[irow];
1009 pi[irow] = -2.0 * weight * value;
1010 }
1011 }
1012 } else {
1013 /* dj but at bound */
1014 if (djval > djFlag) {
1015 statusWork[icol] = 1;
1016 nflagged++;
1017 }
1018 }
1019 }
1020 }
1021 }
1022#ifdef TWO_GOES
1023 }
1024#endif
1025 }
1026#ifdef TWO_GOES
1027 for (i = 0; i < nrows; i++) {
1028 rowsol[i] = 0.5 * (rowsol[i] + rowsol2[i]);
1029 pi[i] = 0.5 * (pi[i] + pi2[i]);
1030 }
1031#endif
1032 if (extraBlock) {
1033 for (i = 0; i < extraBlock; i++) {
1034 double value = solExtra[i];
1035 double djval = costExtra[i];
1036 double djval2, value2;
1037 double element = elemExtra[i];
1038 double theta, a, b, c;
1039 int irow = rowExtra[i];
1040 djval -= element * pi[irow];
1041 /*printf("xxx iter %d extra %d djval %g value %g\n",
1042 iter,irow,djval,value);*/
1043 if (djval > 1.0e-5) {
1044 value2 = -value;
1045 } else {
1046 value2 = (upperExtra[i] - value);
1047 }
1048 djval2 = djval * value2;
1049 if (djval2 < -1.0e-4 && fabs(djval) > djTol) {
1050 nChange++;
1051 a = 0.0;
1052 b = 0.0;
1053 c = 0.0;
1054 djval2 = costExtra[i];
1055 value = rowsol[irow];
1056 c += value * value;
1057 a += element * element;
1058 b += element * value;
1059 a *= weight;
1060 b = b * weight + 0.5 * djval2;
1061 c *= weight;
1062 /* solve */
1063 theta = -b / a;
1064 if (theta > 0.0) {
1065 value2 = CoinMin(theta, upperExtra[i] - solExtra[i]);
1066 } else {
1067 value2 = CoinMax(theta, -solExtra[i]);
1068 }
1069 solExtra[i] += value2;
1070 rowsol[irow] += element * value2;
1071 value = rowsol[irow];
1072 pi[irow] = -2.0 * weight * value;
1073 }
1074 }
1075 }
1076 if ((iter % 10) == 2) {
1077 for (i = DJTEST - 1; i > 0; i--) {
1078 djSave[i] = djSave[i-1];
1079 }
1080 djSave[0] = maxDj;
1081 largestDj = CoinMax(largestDj, maxDj);
1082 smallestDj = CoinMin(smallestDj, maxDj);
1083 for (i = DJTEST - 1; i > 0; i--) {
1084 maxDj += djSave[i];
1085 }
1086 maxDj = maxDj / static_cast<double> (DJTEST);
1087 if (maxDj < djExit && iter > 50) {
1088 //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
1089 break;
1090 }
1091 if (nChange < 100) {
1092 djTol *= 0.5;
1093 }
1094 }
1095 }
1096RETURN:
1097 if (kgood || kbad) {
1098 COIN_DETAIL_PRINT(printf("%g good %g bad\n", kgood, kbad));
1099 }
1100 result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
1101 rowlower, rowupper, lower, upper,
1102 elemnt, row, columnStart, length, extraBlock, rowExtra,
1103 solExtra, elemExtra, upperExtra, useCostExtra,
1104 weight);
1105 result.djAtBeginning = largestDj;
1106 result.djAtEnd = smallestDj;
1107 result.dropThis = saveValue - result.weighted;
1108 if (saveSol) {
1109 if (result.weighted < bestSol) {
1110 bestSol = result.weighted;
1111 CoinMemcpyN(colsol, ncols, saveSol);
1112 } else {
1113 COIN_DETAIL_PRINT(printf("restoring previous - now %g best %g\n",
1114 result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
1115 }
1116 }
1117 if (saveSol) {
1118 if (extraBlock) {
1119 delete [] useCostExtra;
1120 }
1121 CoinMemcpyN(saveSol, ncols, colsol);
1122 delete [] saveSol;
1123 }
1124 for (i = 0; i < nsolve; i++) {
1125 if (i) delete [] allsum[i];
1126 delete [] aX[i];
1127 delete [] aworkX[i];
1128 }
1129 delete [] thetaX;
1130 delete [] djX;
1131 delete [] bX;
1132 delete [] vX;
1133 delete [] aX;
1134 delete [] aworkX;
1135 delete [] allsum;
1136 delete [] cost;
1137#ifdef TWO_GOES
1138 delete [] pi2 ;
1139 delete [] rowsol2 ;
1140#endif
1141 for (i = 0; i < HISTORY + 1; i++) {
1142 delete [] history[i];
1143 }
1144 delete [] statusSave;
1145 /* do original costs objvalue*/
1146 result.objval = 0.0;
1147 for (i = 0; i < ncols; i++) {
1148 result.objval += colsol[i] * origcost[i];
1149 }
1150 if (extraBlock) {
1151 for (i = 0; i < extraBlock; i++) {
1152 int irow = rowExtra[i];
1153 rowupper[irow] = saveExtra[i];
1154 }
1155 delete [] rowExtra;
1156 delete [] solExtra;
1157 delete [] elemExtra;
1158 delete [] upperExtra;
1159 delete [] costExtra;
1160 delete [] saveExtra;
1161 }
1162 result.iteration = iter;
1163 result.objval -= saveOffset;
1164 result.weighted = result.objval + weight * result.sumSquared;
1165 return result;
1166}
1167